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Abstract. Gauge invariance was discovered in the development of classical elcctromagnetism and was 
required when the latter was formulated in terms of the scalar and vector potentials. It is now considered 
to be a fundamental principle of nature, stating that different forms of these potentials yield the same 
physical description: they describe the same electromagnetic field as long as they are related to each other 
by gauge transformations. Gauge invariance can also be included into the quantum description of matter 
interacting with an electromagnetic field by assuming that the wave function transforms under a given 
local unitary transformation. The result of this procedure is a quantum theory describing the coupling 
of electrons, nuclei and photons. Therefore, it is a very important concept: it is used in almost every 
fields of physics and it has been generalized to describe electroweak and strong interactions in the standard 
model of particles. A review of quantum mechanical gauge invariance and general unitary transformations is 
presented for atoms and molecules in interaction with intense short laser pulses, spanning the perturbative 
to highly nonlinear nonperturbative interaction regimes. Various unitary transformations for single spinless 
particle Time Dependent Schrodinger Equations, TDSE, are shown to correspond to different time-dependent 
Hamiltonians and wave functions. Accuracy of approximation methods involved in solutions of TDSE's such 
as perturbation theory and popular numerical methods depend on gauge or representation choices which 
can be more convenient due to faster convergence criteria. We focus on three main representations: length 
and velocity gauges, in addition to the acceleration form which is not a gauge, to describe perturbative and 
nonperturbative radiative interactions. Numerical schemes for solving TDSE's in different representations 
are also discussed. A final brief discussion is presented of these issues for the relativistic Time Dependent 
Dirac Equation, TDDE, for future super-intense laser field problems. 
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1. Introduction 

Advances in current laser technology allow experimentalists to access new laser sources for probing and 
controlling molecular structure, function and dynamics on the natural time scale of atomic (nuclear) motion, 
the femtosecond (1 fs = 10~^^ s) [11^2 and electron motion on attosecond (1 as = 10^^^ s) time scale [3, 4 . 
A new regime of nonlinear nonperturbative laser-matter interaction is leading to new physical phenomena 
such as nonlinear photoelectron spectra called Above Threshold Ionization (ATI) [SJ [BJ [TJ [H] and high order 
harmonic generation (HHG) pillO). In these two new examples of nonperturbative electron response, ionized 
electron trajectories are completely controlled by the electric laser field [3]. The atomic unit of laser intensity 
/o = cEq/Stt = 3.5 X 10^^ W- cm~^ corresponds to the atomic unit (a.u.) of electric field i?o = 5 x 10^ 
V-cm^^, the field at the atomic radius ao = 0.0529nm of the Is hydrogen (H) atom orbit. For an intensity 
/ = 3.5 X 10^"* W-cm~^ = 10~^ a.u. currently used in many high intensity experiments, the electron 
ponderomotive energy at wavelength A = 800nm {lo = 0.057 a.u.), Up = Ij^uP' (a.u.) = 0.77 a.u. = 21 eV, 
exceeds the ionization potential /p = 0.5 a.u. — 13.6 eV of the ground state H atom. The corresponding 
maximum field induced excursion of the electron is a = E/uP — 30.8 a.u. = 1.63 nm, where the Coulomb 
potential becomes negligible and a field dressed electron description becomes appropriate. This new regime 
of nonperturbative radiative interactions has motivated the development of simple but highly predictive 



physical models, such as the recollision model in atoms [TTl [T^l [13] or molecules [5], and the strong field 
approximation (SFA) [3', '5', T| which have become standard models for the advancement and development of 
this new field of nonlinear physics. 

The theoretical description of perturbative and nonperturbative radiative interactions relies on one of the 
most important concepts of physics: gauge invariance. Gauge invariance was discovered in the development 
of classical electromagnetism when the latter was formulated in terms of scalar and vectorial potentials ,14, . 
It is now considered to be a fundamental principle of nature, stating that different forms of potentials yield the 
same physical description, i.e. they describe the same electromagnetic fields as long as they are related to each 
other by gauge transformations |15| . In the quantum description of matter interacting with electromagnetic 
fields, gauge invariance is obtained by transforming wave functions under local unitary transformations, 
resulting in different Hamiltonians in the corresponding TDSE's. However unitary transformations can 
also give rise to "representations" which are not considered as "gauge" transformations (the latter are 
obtained rigorously only from transformations of classical Lagrangians leaving the dynamics invariant [H [15] ) . 
According to the gauge principle (this will be described in more details in latter sections), all physical 
observables are gauge invariant. However, their calculations usually involve the evaluation of gauge dependent 
quantities: for instance, the expression of the electron wave function depends on the gauge chosen. When the 
exact analytical solution of the TDSE is known, it is possible to move from one gauge to the other by a gauge 
transformation (see |16| for instance), and then, all gauges will yield the same physical result. In this case, 
the gauge choice is simply a matter of convenience: it may be easier to solve the TDSE in one specific gauge 
than in others. However, when approximations of gauge dependent quantities are involved (for instance, 
when using numerical methods or perturbation theory) , the gauge invariance of physical observables may be 
lost as the error induced by the approximation scheme may not transform in the same way as the full solution. 
This can be illustrated as follows. Let us consider a single electron in interaction with an electromagnetic 
field expressed in two different gauges such that our system can be described equivalently by the following 
wave functions: 

V-w = ^(i) + i?(i), (1) 

V,(2) =^(2)+i?(2)^ (2) 

where ■0''^'' is the general exact solution of the TDSE, coupled to the electromagnetic field expressed in 
gauge 1, while V'^^^ is given in gauge 2. The two wave functions are related by a gauge transformation G as 
tp^^'' — G^^^\ with the following calculated observables: 

(6) = (6(1)) := (V-WlO^i^W) = (V^'^IO^'^IV^'^) =: (O^^^), (3) 
provided that the observable transform as O'^) — GO'^^^G^^ , which can be verified explicitly in most 
cases. The last equation ([3]) expresses the essence of the gauge principle, i.e. that physical observables 
are independent of the gauge chosen. In Eqs. ([T]) and ©, ■0'^'^'' are approximate wave functions, obtained 
from a solution of the TDSE using approximation methods. Thus, i?(i'2) are the remainder or error of this 
method. Typically, it will be given by i?(i'2) ^ for numerical methods (where h is the grid size and 
n e N*) or i?(i'2) ~ (/" for perturbation theory (where g is a small parameter). If the approximate wave 
functions obey the same gauge transformation as the full solution, that is ^^^^ = G^/'(2) ^ the approximation 
of the observable will also be gauge invariant: 

(6) « (^W|6(l)|V;(l)) = (^(2)|o(2)|^(2)^^ (4) 

However, this is generally not the case because ■(/'(i) and ■0(2) usually transform differently in different gauges 

(6)«(0W|6«iv;«), (5) 

« (0(2)|O(2)|^(2))^ (6) 

but 

(0;(i)|6(i)|0W)^(v;(')|o(')|v;^'^), (7) 

and thus (O^i)) « (0(2)). Therefore, we have lost gauge invariance by approximating the wave function. As 
the approximate wave functions get closer to the exact solution, the observables calculated in two different 
gauges converge towards each other 

(v;(i)|o(i)|0(i)) ^'''"^°> (0(2)|6(2)|^(2)) ^ 
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and of course, gauge invariance is recovered in the limit of the exact solution when R^^'"^") — 0. 

The previous discussion illustrates the fact that approximating a gauge independent quantity (an 
observable) by implementing an approximation of a gauge dependent quantity (the wave function) may 
destroy the gauge independence of the former: the breaking of gauge invariance is an artifact of the 
approximation method. In this paper, we will compare the calculation of many observables in different 
gauges and show that certain gauge choices have better convergence properties towards the exact solution. 
One should always try to choose the gauge which gives the best approximation of the physical quantity under 
consideration. This has been studied extensively using either analytical or numerical approximations. For 
instance, it was remarked by Lamb in his celebrated study of the Hydrogen atom fine structure, that the 
theoretical results obtained in the length and velocity gauges (defined later) differ in perturbation theory [T7] . 
Moreover, the result of the length gauge calculation was in better agreement with experiments, suggesting 
that this gauge would be more "fundamental". This apparent paradox was resolved latter when it was 
observed that the "naive" perturbation theory is not gauge invariant and that special care is required 
to calculate observables in a gauge independent way. Many papers have attempted to clarify this issue 
[IHllTllinilllllSIlllllllIISlIM]- The equivalence of the velocity and length gauge was also demonstrated 
in specific calculations for multiphotons transition probabilities [35] and induced polarization Other 
calculations have shown the equivalence of results in different gauge choices ^71 [25] . Given the status and 
importance of this issue, we will attempt to formulate a consistent gauge invariant perturbation theory using 
general arguments in the following sections. Concerning numerical calculation and gauge independence, it 
is now generally admitted that certain gauge choices are better than others for the solution of the TDSE. 
For instance, it was concluded that the velocity gauge is more appropriate for calculations in dynamical 
laser- matter interactions [IH] ■ Other calculations where gauge choices are compared can be found in [3^ for 
high harmonic generation, HHG. It should be stressed here that these comparisons should be performed with 
great care because the structure of the mathematical equation may change under a gauge transformation 
and the resulting TDSE may require a different numerical scheme. This issue will be discussed in details in 
this work, along with the description of current numerical methods. 

Throughout this work, a single particle system is mainly considered but extension to A^-particle problems 
exists naturally in the literature, and can often be easily deduced. Recent work has examined gauge invariance 
of coupled electron-nucleus dynamics using the time-dependent Hartree-Fock equation (TDHF) in [31j and 
the Time-dependent Kohn-Sham equation (TDKS) density functional theory [32] and also nonadiabatic 
molecular dynamics [33) . However, they are much more involved technically and are outside the scope of our 
discussion. 

It is also assumed that the electromagnetic field is not quantized and treated as a classical field. This 
approximation is valid when the number of photons is large (such as in a macroscopic laser field) such that 
quantum fiuctuations can be neglected [15]. Concluding whether this approximation is true or not for the 
field induced by the particle is a non-trivial task and certainly depends on the dynamic of the system. Also, 
from the practical point of view, this approach simplifies the calculations significantly, which is certainly the 
main reason why it has been so popular among practionners. It is important to note that using a classical 
theory spontaneous emission is not taken into account but must be introduced by "hand" as in HHG. 

This article is separated as follows. In Section [5] the classical and quantum dynamics of particles 
coupled to an electromagnetic field is treated, with a special emphasis on gauge transformations. It will be 
shown how to obtain the TDSE from the classical dynamics of a scalar particle and how it interacts with the 
electromagnetic field. More generally, the gauge principle, which allows to derive the interaction of matter 
with force carriers will be described. In Section[3j the dipole approximation is detailed along with the regime 
where it can be used. The gauges commonly used in laser-matter interaction are reported in Section|4l More 
precisely, we will consider the length, the velocity and the acceleration gauges. The transformations allowing 
to change the gauge representation are also described. SectionOcontains details on analytical approximations 
used to evaluate quantum transition amplitudes. For instance, a general description of perturbation theory 
and the strong field approximation, with respect to gauge transformation, is included. The gauge invariance 
of these techniques is also shown in some specific examples. Section [6] contains a discussion of numerical 
methods for the solution of TDSE's. It is emphasized that the mathematical structure of the TDSE depends 
on the gauge and thus, the numerical scheme used to solve the equation should be chosen by taking this 
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fact into account. Finally, gauge transformations in relativistic quantum mechanics, which is relevant for 
ultra-intense laser field, is considered in Section [T] 



2. Classical and quantum dynamics of particles coupled to an electromagnetic field 

In this section, we recall Maxwell's equations and the notion of gauge transformations and gauge choice, as 
well as some basic informations about the Hamiltonian and Lagrangian formulations of classical mechanics. 
The latter will then be applied to the case of a single particle interacting with a classical electromagnetic field. 
The corresponding quantum dynamics will be derived at the end of this section using the usual canonical 
quantization scheme. This allows to obtain a TDSE describing the quantum dynamics of a particle coupled 
to an electromagnetic field: this equation is the basis of this paper. Finally, we will discuss the derivation 
of the TDSE from the point of view of the Gauge principle and show how a general symmetry principle can 
be used to obtain an interaction of matter with an electromagnetic field. This non-exhaustive survey will 
allow us to set the framework. This first part is a summary of some key sections of |15' (in particular we 
use similar notations). 



2.1. Maxwell's equations 

The electromagnetic field can be characterized by two physical quantities: the electric field E(a;) and the 
magnetic field B(a;), where x := (x, is a four-vector. These are vector fields which exist only in the 
presence of electric charges: an electric field is produced by a stationary charge while the magnetic field 
is induced by charges in motion (currents). Their dynamics and the relation between them is governed by 
Maxwell's equations. Denoting the electromagnetic field by £ :— (E, B), microscopic Maxwell's equations 
in non- Gaussian units are written as: 

1 

-p(x,i). 



V •E(x,i) 

V •B(x,i) 

V X E(x,i) 

V X B(x,i) 



0, 

-^B(x,t), 
1 d 



(9) 



= 37^E(x,i) + — j(x,t). 



'■dt 



where p is the charge density, j is the current density, Eq is the vacuum permittivity and c is the velocity of 
light. When we are considering a system of I point-like particles of charges qi located at ri(t) at time <, the 
charge density is given by [T5] 



p(r,t) =^(Z,^(r-r,(t)), 



(10) 



and the current density by 



j(r,0 - Vg»r,(t)5(r-r,(i)). 



1=1 



(11) 



Of course, there exists a generalization of these last two formulas to the case of a continuous distribution 
of charges |14j . but it is not required in this paper. The particle positions are determined from the 
Newton-Lorentz equations, which describes the non-relativistic dynamics of point-particles immersed in 
an electromagnetic field. Thus, the classical particles trajectories are solution of: 

d2r,(t) 



(12) 



- g.[E(r,(i),i)+v,(i) X B(r,(t),t)], 

ri(to) = R-i, 
.v,(to) = V„ 

for i = 1 • • • ^, and where we consider the electromagnetic field as external. Here, Vj are the particle velocities, 
is the initial time and Ri, are the initial particle positions and velocities, respectively. 
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From the divergence free equation V • B = implying the non-existence of magnetic monopoles, and 
the Maxwell-Faraday equation, we deduce the existence of vector and scalar functions A and U such that 

B = V X A, 
d 



E 



dt 



(13) 



where A is the vector electric potential and U the scalar electric potential. Maxwell's equations can be 
rewritten in terms of A and ?7 as a second order (in time and space) wave equation: 



AC/ + - V • A = p. 

Ot Eq 



1 d_ 



U 



(14) 
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This system of equations is equivalent in principle to Maxwell's equations: it should give the same 
electromagnetic field £. However, we will see in the following that this is not quite the case because the 
potentials are not defined uniquely and thus, we need another condition to obtain the appropriate dynamics. 
This can be understood in the following way. Let us consider a regular scalar function F{r,t). It can be 
easily shown from that the electromagnetic field £ is unchanged (using the fact that V x VF = 0) by 
the so-called gauge transformation: 

A ^ A' = A + VF, 



U 



U' 



d 



(15) 



As a consequence £, is not uniquely defined by the potentials: there exists an infinite number of potentials 
related by gauge transformations that yield the same electromagnetic field. An additional condition, called 
a gauge condition, allows to fix these superfluous "degrees of freedom" and to determine a unique definition 
of potentials; this procedure is called gauge fixing. There exists in principle an infinite number of these 
conditions but from a practical point of view, only a certain number of gauge choices are generally utilized 
because the resulting equations are simpler or have certain interesting properties. In the framework of 
electrodynamics coupled to non-relativistic classical particles, the most popular choices are the Lorentz and 
Coulomb gauges: 



• Lorentz Gauge: 

Consists of imposing the condition 

d 

— [/ + c^V • A = 0, 
ot 

and Maxwell's equations become 



c 

—P, 

£q 
1 

— j- 

£0 



(16) 



(17) 



This gauge has the merit of being manifestly covariant, which is a very useful property when one is 
interested in symmetries of Maxwell's equations or in the covariant perturbation theory. It is to be 
noted that scalar U, and vector A potentials are decoupled in this gauge, as well as charge and current. 

Coulomb Gauge or minimal coupling: 

Consists of imposing the condition V • A = 0, and as consequence Maxwell's equations are rewritten 
AU = --P, 



£0 



1 



d 



(18) 



c'AA = -j-V-C/. 

So ot 



I It is interesting here to note that it took almost a century to deduce this condition |34| . 
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There exist many other possibilities, such as the light-cone gauge, temporal gauge, axial gauge, Fock- 
Schwinger gauge and Poincare gauge (for an exhaustive enumeration and definitions, see [15i 134] ) . Also, 
it is interesting to note that the null divergence of the vector potential A has interesting connections to the 
incompressibilty condition in fluid dynamics. This is discussed in more details in Section [51 



2.2. Lagrangian and Hamiltonian formulation 

This section is not an exhaustive presentation of the notion of Lagrangian and Hamiltonian operators. 
However simple facts are important to recall, in particular the least action principle, which lies as the basis 
of classical and quantum mechanics. Two cases will be treated: the discrete and the continuum cases. The 
former is used for the description of point-like particles. The latter is important when one is interested in 
the dynamics of a field variable such as the electromagnetic field, the velocity field in fluid dynamics and 
even the quantum wave function. 



2.2.1. Discrete case: particle-like systems We denote by (ri)i=i,... .f the set of trajectories of £ particles 
of mass (mi)i=i^... J. In the Lagrangian formalism the search of trajectories is equivalent to solving an 
extremum problem for the action S, between times ti and t2 (see [351 136] ): 

S:^ I L{vi{t),--- ,re{t),vi{t),--- ,re{t),t)dt, (19) 

V(ri(t),... ,r,(i)) -y(ri(t),... ,r,{t))dt, (20) 

where L := K — V is called the Lagrangian, with K = "^i(ri)^/2 the kinetic energy and V the potential 
energy. Above and in the following, the notation d denotes a time derivative of a. It should be noted here 
that the Lagrangian depends only on the variables r^, their time derivative r (velocities) and possibly on 
time, but not on the acceleration r, which is not included. The variables ri,ri are the dynamical variables 
and they completely specify the state of a classical system (this is the reason why the Lagrangian does not 
depend on higher time derivatives) . The explicit time dependence of the Lagrangian is included to describe 
external forces acting on the dynamical system under consideration. In this latter case, it can be shown that 
the energy is not conserved. 

The equations of motion are then obtained from the least action principle, which states that the particle 
paths minimize the action, that is: 6S — J^^ 6L — where 6L is the functional differential: 

^^-Efe^r. + -.r.). (21) 

2—1 

Assuming that the coordinate variation vanishes at t = ii,2, the Euler-Lagrange equations can be obtained 
[35J, for aU i = I,-- • ,£: 

dL d^dL\ 



dvi dt\dvJ 



0. (22) 



These equations allow us to obtain a description local in time (equation of motion) from a global in time 
principle (least action principle). 

Another important quantity can be obtained from the Lagrangian: the conjugate momentum. It is 
given, for the ith particle, by 

dL 

P. = ^. (23) 

It should also be noted that the conjugate momentum will be especially important when the theory is 
quantized and it will happen that Pi is not always equal to mr^. 

The last equation suggests that it is possible to obtain a description of the dynamics in terms of momenta 
p and coordinates r, instead of the velocity r used in the Lagrangian formulation. This is the Hamiltonian 
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formulation, which is related to the Lagrangian by a Legendre transformation: 

Hin,- - ,re,pi,--- ,pi)=Y,^iPi-L = K + V. (24) 

i=l 

Thus, the Hamiltonian represents the total energy (kinetic + potential energy) of the system. From the 
above equation it follows that if the Lagrangian is time-independent, the Hamiltonian or the total energy is 
conserved and thus, is a constant of motion. 



2.2.2. Continuous case: classical field theory In this section, we recall some basic facts about Hamiltonian 
and Lagrangian operators in classical field theory. The main diflFerence with the preceding section is the 
fact that now, the object under study are fields, that is quantities which take a value at each point of space 
time. The latter come in various form depending on their transformation properties: scalar fields, vector 
fields or tensor fields for example. Also, they can be used to describe many physical entities such as the 
electromagnetic field, flow velocity in fluid mechanics, temperature distribution in a material, etc. Their 
dynamics can also be formulated in terms of Lagrangian and Hamiltonian mechanics, which is the subject 
of this section. 

First, let us consider the dynamical variables given by (l)a{x) and di(j)a{x), that is the field under 
consideration. Here the index a= 1, • • • , n is an integer which denotes one of the n field vectorial components, 
that is: 



(t)i{x) 



(25) 



_<t>n{x)_ 

Also, the argument is x = (t,x) such that (j) has a value over x G and t € [ti,tf]. Finally, the index i 
denotes a derivative with respect to a space-time coordinate (when considering the Minkowski metric or in 
other words, when looking at relativistic systems, the usual notation is to have z = /x where /x is a Lorentz 
index). 

Now, since the Lagrangian is a function of the dynamical variables, we can write the field action as 



= / 

Jti JR^ 



<fxjC[^a{x),di(l)a{x),x], 



(26) 



where £ is a Lagrangian density. This action is a scalar and is a generalization to the continuous case of 
the action for the discrete case. It can actually be derived from the discrete case in the limit of an infinite 
number of degrees of freedom and by assuming local interactions. The least action principle is unchanged 
in this procedure and can still be used to compute the equations of motion. The latter states that the field 
minimizes the action such that 6S = S[(j)'] — <S'[(/)] =0 where (l>'a{x) = (j>a{x) + S<j)a{x). In other words, the 
action is stationnary under the perturbation 6(pa{x). So under this infinitesimal variation, the action changes 
according to 



6S 



dt 



x 



dC{x) 



di 



dC 



d(l)a{x). 



d(j)a{x) d[di(t)a{x)] 

Of course, by requiring stationarity of the action, we get the usual Euler-Lagrange equation of motion: 
dC{x) dC{x) 



= 0. 



(27) 



(28) 



d(j}a{x) ^'d[di(pa{x) 

Here, Einstein notation convention is assunicid, that is repeated indices are summed. 

The conjugate momenta and the Hamiltonian density can also be defined in a similar way as in the 
discrete case. They are given by 

dC 

"a (a:) = TT-T^, (29) 



d(t)a{x) 
^^a{x)<j)a{x) 



c. 



(30) 
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Finally, it is often convenient to use the Lagrangian of a contimious system obtained from the Lagrangian 
density as L = J^3 (Px£{x). This notation will be used frequently in the following sections. 

2.2.3. Symmetry transformations In both the discrete and the continuous cases, the physical system may 
have symmetries, that is a set of transformation which leave the dynamics invariant. In the Lagrangian 
formulation, it can be shown that under such symmetry transformation, the Lagrangian is unchanged, up 
to a divergence term. Symmetries are very important in all areas of physics because they are related to 
conserved quantities via Noether's theorem. 

First, we will look at the possible transformations that can be implemented on the Lagrangian. These 
exist in two varieties (in the discrete case, only the first one can be implemented): 

(i) Transformations on coordinates: 

Examples of these are the Lorentz and Galilean transformations, which include translations and 
rotations. 

(ii) Transformation on the fields: 

Examples of these are the phase transformation of the wave function {ip — > tj/ = e^^ijj ). Note that in 
the following, we will consider a set of invertible transformations which depend only on the field itself 
(not its derivative). 

Mathematically, a general way of writing these two transformations is to define a linear mapping Ta as: 

T ' ^ '-^A^ (31) 

^ ' (l)a{x) (pA.ai^A)^ 

where A is a continuous parameter (which may depend on spacetime) such that when A = 0, the 
transformation is the identity. Then, it can be shown that a condition to confirm that Ta is a symmetry 
transformation (for the continuous case) is that 

J{x,XA)jC.[4>A{XA),d^(l)A{XA),XA] = jC[(f>{x), di(f>{x), x] + diF^^[(f){x),x], (32) 

where J{x,xa) is the Jacobian of the coordinate transformation. In the discrete case, one finds that this 
symmetry condition is written as 

L = L+jJ{ri,t), (33) 

where /, F arc arbitrary functions. If these conditions are fulfilled, the transformed and the initial Lagrangian 
will lead to the same dynamics (the equation of motion will have the same form). 

2.3. Lagrangian of the electromagnetic field 

In this section, we will look at the Lagrangian density describing the electromagnetic field. It is convenient 
for this discussion to use the manifestly covariant formulation of electrodynamics. In this case, we define the 
four vector potential and current as 

A^^{x):^iU(x)/c,Aix)), (34) 
J^{x) := {cp{x),3{x)). (35) 

These quantities arc contravariant tensors which are related to their covariant counterparts by = r]^'^A^, 
where rj'^'^ = diag[l, — 1, — 1, — 1] is the Minkowski metric. It is also convenient to define the antisymmetric 
field strength tensor as 

F^"'{x) := d^'A-'ix) - d^A'^ix). (36) 

Then, within this formulation, Maxwell's equations can be written in a very compact and manifestly covariant 
(invariant under Lorentz transformations) form: 

d^F^'-'ix) = -^r{x). (37) 
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Also, the conservation of current is written as 

c)^J^(a;)=0, (38) 
and finally, the gauge transformation is 

A^'ix) A'^'ix) = A^'ix) + d^'Fix). (39) 

We would like to stress that the covariant formulation is equivalent to the set of Maxwell's equations presented 
before. However, it is independent of the choice of referential frame: the latter are related by Lorentz (or 
Poincare transformations) . 

From these equations, it is possible to obtain the corresponding Lagrangian density for the 
electromagnetism sector. It is given by 

£E&M(A^VA^,x) :=/:kin(^^VA^x)+/:i„t(v4^,V^^a;), (40) 



■F^,{x)F^''ix)~J'^ix)A^{x). 



(41) 



It is a straightforward calculation to show that the Euler-Lagrange equation for this Lagrangian density is 
given by Maxwell's equations. 

The coupling of the electromagnetic field to matter is done through the source terms J^A^ which contain 
the charge density and current. In a classical setting, this would be related to the charge position and thus, 
another term should be included in the Lagrangian to take care of the particles dynamics. This is done in 
section ^TM In the quantum setting however, the tensor is related to the wave function of the particle 
under study (in our case, a single electron). This is described is Section ^IM 



2.4- Lagrangian of a particle in an electromagnetic field 

The Lagrangian for a system of £ non-relativistic free particles is given by Lf-part l/2X]i=i '^i'^'i- To 
introduce its coupling with an electromagnetic field, we can add the electromagnetic field Lagrangian Le&m 
and we get 

Llp+S '■= -^£-part + -^E&M, (42) 



1 ^ t 



(43) 



i=l 



' i=l 



{-VUix) - k{x)f - c2(V X A(a;))^ 



i(x)-A{x)-p{x)U{x) 



(44) 



S{x)-A{x)~p{x)U{x) , (45) 



where the current j and density p{x) in the interaction terms of £e&m are given by Eq. ([Tl]) and (jlOp . 
respectively. The Euler-Lagrange equation of motion are given by Maxwell's equations and by the Newton- 
Lorentz equation as the dynamical variables are r, r. A, C/. 

We now discuss the effect of a gauge transformation via F (IT5|) . on this Lagrangian where we set ^ = 1, 
that is for the single particle case. Replacing (A, [/) by (A',U'), leads to a new Lagrangian L: 



L ^ L:=L + j\v ■{iF) + ^^{pF) 



Due to charge conservation 
dp 



(46) 



(47) 
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and the divergence theorem (/V • {jF)d^r ~ 0), we easily deduce that 



(48) 



L+— pFd?v. 
at J 

The new L and old Lagrangians L are then equivalent as they obey the symmetry condition in Eq. psp . 
This transformation is called a gauge transformation of the first kind, according to Pauli ^37j . The previous 
Lagrangian was written in a general way such that gauge invariance is explicitly satisfied. However, it 
contains redundant degrees of freedom and as discussed earlier, this can be discarded by the gauge fixing 
procedure. Thus, it is possible to write Lagrangians for specific gauge choices. For instance, in the Coulomb 
gauge (B = V X A with V • A = 0), the Lagrangian (j45|) can be written: 



Vc 



- c2(V X Af +j • A 



£0 

2 ' ' 2 

and the corresponding Hamiltonian, is then given by (using (|24p) 



-c2(V X A) 



(49) 



(50) 



There are a few interesting remarks to make about this Lagrangian: 

• It is not gauge invariant, rather, it is obtained from the gauge invariant Lagrangian by gauge fixing. 
Thus, it is valid only for the Coulomb gauge choice. 

• The electromagnetic field dynamics does not depend on the scalar potential U. The gauge fixing 
procedure (V • A = 0) allowed to eliminate this degree of freedom. 

• The Coulomb potential Vc of the particle appears naturally from the term involving the charge density 
[15| . Thus, in this gauge, the field of the particle is simply given by the usual Coulomb law. 

• It is not manifestly covariant because the gauge condition is not invariant under Lorentz transformation. 

For these three reasons, the Coulomb gauge is a very popular choice to describe laser-matter interactions. 

We would like to conclude this section by considering the special case of a single particle (£ — 1) subject 
to an external electromagnetic field where we assume that the field is not part of the dynamical system. 
This approximation is often used to simplify the calculations. The simplified Lagrangian becomes 



^p+cxt 



d-'r 



j{x) ■ Ae{x) - p{x)Ue{x) 



(51) 



where Ue,Ae represent the potential of an external electromagnetic field £e '■— (Eg, Be). In this model, 
the Euler-Lagrange equation are given by the Newton-Lorentz equation for the particle and there is no 
backreaction of the particle on the field. The conjugate momentum is p = rnv + qA,, while the Hamiltonian 
is 

1 



p+ext 



2m 



[p-qAe] +qUe 



(52) 



This Hamiltonian is then used in the special case where it is assumed that the backreaction on the 
electromagnetic field is negligible. Throughout this paper, we will refer to this case as the external field 
approximation (minimal gauge). 



2.5. Quantization 

The equations of motion obtained in the last section can be quantized in the usual canonical quantization 
where essentially the position and conjugate momentum become operators. This method attempts to 
quantize a classical system while keeping its main properties such as symmetries. We would like to stress 
again that in this work, we are not quantizing the electromagnetic field: only the single particle described 
by the Newton-Lorentz equation will be quantized. The canonical quantization states that the classical 
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dynamical variables, that is the position r and conjugate momentum p, becomes operators with the following 
commutation relations: 

[f„f,]-0, (53) 
[p^,Po] =0, (54) 
= ihSij. (55) 

Throughout this work, we will work in the "position representation" where |r) is a vector in the Hilbert 
space describing the quantum state of our system. The position operator has the property that f |r) = r|r). 
From this result and the commutation relation, it is straightforward to obtain that p = —ihV (in the free 
case). 

A general state is obtained by the linear superposition ~ J (Prip{r)\r). The wave function is a 
projection of such a state on the position state, that is ip{r) (r|V')- The dynamics of the wave function is 
given by the TDSE: 

ih^m)) = mm)). m 

where H is the Hamiltonian operator for the system under consideration. Projecting this equation on position 
space, we get the usual TDSE in coordinate space: 

d 

in-i;{t,r) ^ H{t,r)i^{t,r), (57) 

where the quantum Hamiltonian is obtained from the classical Hamiltonian by using the following 
prescription: p — ^ —ihS/. This procedure will be used in the rest of this paper to obtain wave equations 
describing the quantum dynamics of a single particle interacting with different electromagnetic fields. 



2. 6. Time dependent Schrodinger equation coupled to an electromagnetic field and the gauge principle 

In this section, we consider the gauge invariance of the TDSE coupled to an electromagnetic field from the 
Lagrangian viewpoint and the gauge principle. It should be noted here that just like the electromagnetic 
field, the wave function is also a field (a complex scalar field) in the sense that it is a function that is 
defined over all M.^. Therefore, its dynamics can be understood from a variational principle analogous to the 
treatment of the electromagnetic field. This will be presented in the first section. Finally, the coupling of 
the two along with gauge invariance and physical consequences will be proven. 



2.6.1. Schrodinger scalar field The free Schrodinger equation for a particle of mass m is given by 

ihdtip(x) = -h^—t/jix). (58) 
2m 

Then, it can be shown that the free Schrodinger equation is given by the Euler-Lagrange equation of the 
following Lagrangian density: 

Cs{^,di^,r,dr,x) = ihdt\i>ix)f - —\^^{x)\\ (59) 

Zm 

It should be noted that since the field ip is complex, both i/; and ip* should be considered as dynamical 
variables [15j . leading to two different Euler-Lagrange equations which are complex conjugate of each 
other. Also, the Schrodinger Lagrangian is not invariant under Lorentz transformations because of course, 
it describes a non-relativistic particle. It is however invariant under the Galilean transformations which 
relates different referential frames in the non-relativistic setting. There exists relativistic generalizations of 
the Schrodinger wave equation, such as the Klein-Gordon (spin-0) and the Dirac (spin-1/2) equations, which 
are invariant under Lorentz transformations. 
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2.6.2. Minimal coupling prescription and gauge invariance Now, we would like to couple the "matter field" 
described by the Schrodinger Lagrangian with the electromagnetic field. This can be achieved by imposing 
a symmetry on the Lagrangian £5. Let us assume that the theory describing our system is invariant under 
local phase transformations {U{1) symmetry), that is: 

i^{x) i}'{x) = e*^(^)/'"V(a:), (60) 

such that C'g = £3. However, an explicit calculation shows that Cg ^ £5. Rather, we have 

C's{i^,d^,r,dr,x) = thdtmx)\^ - ^UnW - VF(x))VXx)p + dtF{x). (61) 

The only way to cancel the extra terms in Eq. (j6ip is to add a new field with appropriate transformations. 
This new field is the electromagnetic field and it is added via the minimal coupling prescription, that is 
partial derivatives are replaced by: 

^ d" + eA''{x). (62) 

This is complemented by adding the kinetic term F^'^ F^^i, such that we obtain a dynamical field. Then, the 
Lagrangian for the coupled system becomes 

-Cms :=£Ms(V',aV',^^9V^^^VA^x), (63) 

= [ihdt-eU{x)]\tP{x)\^-^\[-ihW-eA{x)]i^{x)\^-^-^F^^x)F^,{x). (64) 

2m 4 

This represents physically the dynamics of the matter field with its "own" electromagnetic field, that is the 

magnetic field generated by the particle itself. However, it is also possible to add an external electromagnetic 

field by letting A^^ A^^ + A^^j . The corresponding Euler-Lagrange equation then becomes 

ihdtip{x) = ^ [-iftV - e{A{x) + Ac^t{x))f ip{x) + e[U{x) + Ue^tix)]4>{x), ^^^-^ 



d^Ft'^x) ^ J-'ix), 



where 



e|i/'(a;)|2 for = 0, 



[[-ihW -e{A{x)+ A,^t{x))]\^{x)\^ ior pi = 1,2,3. 

Before continuing further, we summarize the above derivations. We started with a Lagrangian describing 
matter: the Schrodinger Lagrangian in Eq. (|59|). Then, we imposed a local C/(l) symmetry (note that U{1) 
symmetry with a space-independent phase parameter is related to charge conservation) to this Lagrangian. 
The consequence of this is that we had to add a new field, the electromagnetic gauge field, such that the 
symmetry is obeyed. This field obeys the gauge transformation properties. This whole procedure is an 
example of the gauge principle which is used in many field of physics to obtain interaction terms between 
matter and force carriers. For example, the theory of Strong and Weak nuclear interaction are based upon this 
very principle, albeit on its non-abelian generalization (the symmetry groups are local SU{S) and SU{2) in 
these cases). Therefore, it is generally believed that gauge symmetries are one of the fundamental organizing 
principles of nature. 

To summarize, the gauge transformation T\ can be written as 

. ^{x) ^ V'(a^) = e*^(-)/V(x), ..-^ 
Af'ix) A'^ix) = A^(x) +a^F(a;). ^ ' 

where A is an arbitrary function. This was shown to be a symmetry transformation because it obeys Eq. 
(El. 



2. 7. Gauge invariance, unitary transformation and quantum observables 

The goal of this section is to show that although gauge transformations are unitary, not all unitary 
transformations correspond to a gauge transformation. This may occur when the transformation parameter 
F depends on dynamical variables such as r, etc. Also, the fact that transition matrix elements are invariant 
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under any unitary transformations is discussed. This is crucial for our analysis because it states that physical 
observables are the same, for any unitary transformations of the wave function. The corollary to this is that 
observables are invariant under gauge transformations. 

To describe these results, we will start by looking at the invariance of quantum observables under unitary 
transformations. Let us consider a transformation operator defined by 

U{t) := eyiY>{iF{t)/h), (68) 

where F can be any space-time dependent operator. Of course, this transformation is unitary: |Z//(t)p — 1. 
The operator lA acts on the wave function as IV'^^HO) = ^(OIV-^^HO), while the wave functions obey the 
following TDSE's: 

It can be verified by using these TDSE's, along with the definition of unitary transformations, that the 
Hamiltonians in the two representations are related by 

= UH^^)u^ + ih^K (70) 
dt 

Note here that this is derived by assuming that [dtU,U] ^ 0. The average energy is then given by 

dt 

Therefore, the average energy is not invariant under a general unitary transformation if the last term of the 
last equation is non-zero. However, it can be easily deduced that physical observables, which are given in 
terms of transition amplitudes, are invariant under unitary transformations, even if the Hamiltonian does 
not. This result can be obtained as follows. 

Denoting U'^^^^'>{t,tQ) := f exp(-i //^ H^^''^\t')dt' /K) the evolution operator for ift^'^) (where f stands 
for the time-ordering operator) see for instance |38], such that for t ^ to, 

\^^'\t))^m'\t,to)\iJ^'\to)). (72) 

Then, we have 

\,p^'\t))^um^'\t)) 

=u{t)m'\t,t,M^'Hto)} 

^U{t)U^'\t,to)uHtoM^^\to)) 

= U('\t,to)\iJ^^\to)). (73) 

Thus, we can deduce that the evolution operator transforms as 

U^'\t,to)^Um^'Ht,to)uHto) (74) 

under a unitary transformation. As a consequence, the transition amplitudes from |'(/;^^^(to)) to lip^^^t)) 
obey: 

{i/'\t)\U^'\t,to)\i/'\to)) = {r^'^'Ht)\U^'Ht,toM^'Hto))- (75) 

This equation states that the transition amplitudes are equal and are invariant under a unitary 
transformation. This is a very important result because most physical observables can be obtained from the 
transition amplitudes and consequently, these observables are also invariant under unitary transformations. 
For instance, in HHG experiments, one can measure the amplitude and phase of each harmonic electric field, 
and these are related to photon emission transition moments, which can result in direct tomography of wave 
functions [^1155]. 

We can now specialize the general unitary transformations to the special case of gauge transformations. 
Then, the unitary transformation takes the form: 

G(t) =exp(i|F(r,A,[/,i)). (76) 
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In this case, F and ^ commute and the transformation of the Hamiltonian is given by 

H(2)=GHWGt-^. (77) 

It was shown in previous sections that under this transformation, the Lagrangian L*^^-* of a particle 
interacting with electromagnetic radiation transforms into an equivalent Lagrangian via 

L(2)(X,X) = L(i)(X,X) + ^qF(r, A,t/,i), (78) 

where we set X — {r,A,U). This is the more general unitary transformation that yields a Lagrangian 
respecting the symmetry condition in Eq. Note that the function F is independent of r, A or U, to 

avoid any dependence in r, A or U in the Lagrangian as these quantities are not required to specify the 
dynamics of the system. The consequence of adding these terms is that the new Lagrangian L*^^^ would not 
describe a physical system obeying Hamiltonian mechanics. Thus, certain unitary transformations change 
the classical dynamics of the system, such as 

U(t) = exp («|i^'(r, r, A, U, A, [>, t)) , (79) 

although they leave transition amplitudes invariant. These are not gauge tranforniations because they may 
change the form of the Schrodinger equation. An example of these are the Kramers- Henneberger (also called 
Bloch-Nordsieck [3D]) transformations, presented below, which allow to obtain the acceleration gauge from 
the length gauge. 

3. Long wavelength or dipole approximation 

In this section, we describe the long wavelength approximation which is relevant when the wavelength of 
the electromagnetic field A is much larger than the dimension of the system Ac, that is A ^ Ac. In that 
case, the spatial variation of the electromagnetic field over the size of the quantum system is very small and 
thus, it can be neglected. This approximation is often used in laser-matter interaction as it simplifies the 
calculations significantly. 

More precisely, we consider a quantum system centered in r = and an electric field with a space-time 
dependence given by E = E(a;i — kz). Here, we assume that the z-axis of the coordinate system is in the 
direction of the wave propagation and thus, we define the wave number k := |k|. This form of the electric 
field is relevant as E is a solution of a wave equation, which has plane wave solutions with this space-time 
dependence (for instance E = Eq cos(a;i — kx)). A general solution can be written as the linear combination 
of plane waves. Then, making a Taylor expansion as in |41j . 

d 

d{ujt — kz) 

This last equation can be re-written in two different but equivalent ways: 



E(cj< - kz) - E(wt) - fcz^^-^ — E(a;i - A:z)|^^g . (80) 



d 

E(ojt - kz) - E(ojt) + z— E(ojt - kz)\ „ , (81) 
dz 

E(ujt - kz) - E(w<) - — -^E(w<). (82) 
uj ot 

The first one allows us to understand the limit of the long wavelength approximation. Indeed, from this 
equation, we obtain the change in the electric field over the size of the system: AE ~ \cdzE(ujt — kz)\z=Q ^ 
Acfc|E| = Acw|E|/c. To neglect the space variation of the electric field, we need the condition AE ^ |E| and 
thus, we obtain: 

AcCJ < c. (83) 
When the laser frequency obeys this condition, the long wavelength approximation can be used. 
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The second equation allows us to get the vector potential in the velocity gauge associated with the 
electric field, in this long wavelength approximation. It is given by 

A{ujt - kz) - A{ijjt) + -E(wt). (84) 
c 

This form of vector potential with the second term neglected will be used extensively in the next sections. 
We can focus now on the magnetic field, which is given, within the long wavelength approximation, by 

B - iv X [zE(wt)]. (85) 
c 



This is obtained by using the expression of the vector potential in Eq. (184)) . Thus, to first order in the 
approximation (if A — A{ujt)), we have B = 0. However, even if the condition in Eq. (j83| is fulfilled, if the 
electric field is strong enough, the magnetic field cannot be neglected by virtue of Eq. (|85p . Therefore, we 
need another condition for the validity of this approximation. It can be obtained from considerations using 
classical mechanics. The condition is that the electron displacement S due to the magnetic force should be 
smaller than the system size, that is 5 <C Ac. From the Lorentz equation, we know that the magnetic force 
is |F,„ag| ~ e|v||B| = m|a,„ag|, with amag the acceleration due to the magnetic field. The velocity of the 
electron can be estimated from the electric force as |Fciec| = e|E| = m\si\. The typical time for the electron 
acceleration is one cycle: 6t ~ uj~^, so we get that |v| ~ e|E|/(ma;). Then, using the fact that |B| ~ |E|/c 
and |aniag|'5i^, we get the condition 

E« J^, 86 

e V A 

where A = 2t:c/lu is the wavelength of the electromagnetic radiation. This condition on the electric field is 
required to neglect the effect of the magnetic field so that we can use the long wavelength approximation to 
first order. The latter is usually referred to as the dipole approximation: higher order terms corresponds to 
a multipole expansion. 

4. Length, velocity and acceleration gauges 

In this section we detail the main gauge choices commonly used in non-relativistic laser-molecule interaction. 
The starting point of this discussion is the Hamiltonian in the external field approximation written in the 
Coulomb gauge, and given by 



^(Coul) ^ J_ 

2to 



2 



p-qA + Vc, (87) 



where p — —ihV. As discussed previously, this Hamiltonian is obtained from the general gauge invariant 
Hamiltonian by imposing the gauge condition V • A = 0. 

J^.l. Dipole approximation 

To obtain the Hamiltonian in the velocity gauge, we use the dipole approximation described in the last 
section. This allows us to neglect the space dependence of the vector potential and we get 

H^""^ = ^(p-gA(")(c.i))' + K(r), (88) 

which is the Hamiltonian in the velocity gauge [§|. 

Implementing the following unitary transformation 

G'"') {t) := exp {iqA{ujt) • r/^) , (89) 

along with the corresponding gauge transformations of potentials, lead to the length gauge representation: 

= — p2+qr-E(wi) + K(r). (90) 
2m 

§ It should be noted here that the definition of velocity gauge used in laser-matter interaction is different from the one found 
in other contexts, such as in T6 . In the latter, the velocity gauge is a gauge condition which interpolates between the Coulomb 
and Lorentz gauge condition. 
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Still working under the dipole approximation and introducing the following unitary Kramers- 
Henneberger's transformation {AT, also called Bloch-Nordsiek, just introduced for the Dirac equation 
00111113111]: 

G('°'(0 :=exp(- -^gi- /" A(0, s)ds), (91) 



the length gauge Hamiltonian is transformed into 

Hi-} = l-p'^ + ^A^{0,t) + vJr+^ fA{0,s)ds). (92) 
Zm Zm \ m J / 

Due to the presence of r, the unitary transformation does not correspond to a gauge transformations, although 
it preserves the value of observables at all times. Physically, it corresponds to a change of incrtial frame 
as the operator G^'"-* is a translation operator. Although this choice is often referred to as the acceleration 
gauge, it is actually not a gauge choice. For this reason, it has been called the acceleration frame by some 
authors [21], which is certainly a more precise terminology. The main interest of this transformation is to 
remove the transport operator r • E from the Hamiltonian: the inertial frame "follows" the classical motion 
of the particle (without the Coulomb field). The counter-part is that the Coulomb potential is moving jl^, 
which requires a special treatment in the numerical calculations. Its main advantage is that the radiative 
interaction becomes negligible at large distances, such as in ionization, whereas the length gauge with action 
r • E diverges for |r| — > oo. 

4-. 2. First order approximation 

|kT|2 

In the first order approximation, we neglect terms of higher order in 7^ — and keep the first correction 

9 

to the dipole approximation. This may be useful when we are interested in the effect of magnetic field and 
thus, this approach allows for a beyond-dipole approximation study. Starting from the Coulomb gauge and 
keeping the second order term in the long wavelength approximation, we get 

^(2.") = J-(p-gA(")(cji))Vik-r(p-gA(")(wO) • E(wi) + y,(r), (93) 

|k • 

where we neglected terms of higher order in — . Note that the term A(ujt) ■ Eicut) obtained in this way 

is a drift term induced by the magnetic fields. 

The length gauge and acceleration frame representation are obtained in the same way as in the dipole 
approximation. The corresponding Hamiltonians are given by 

^1 i 

jj{2,i) ^ ^2 + _ _(k . r)V) • E,(Lut) + VJr) 

2m ^ c ' 

^(2,a) ^ ^ ^j^2 (^^) ^vJv+^f A(0, S)ds) , (94) 

2m 2m \ m J / 

- r • E(r, t)+qr- VC/(0, t) + — [ A(0, s)ds ■ WU{0, t). (95) 

mj 

4-. 3. Beyong long wavelength approximations 

When the fields are sufficiently strong or when the frequency is large enough, the long wavelength 
approximation cannot be used. However, as shown in [46], it is still possible to derive equivalent Hamiltonians 
in the length gauge and acceleration frame, although their form is more intricate. Starting from the 
Hamiltonian in the Coulomb gauge, we can perform the unitary transformation given by 

Q{^\\,vi) _ ieA{7^) ■ r) , (96) 
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where we defined rj := — k • r. From this unitary transformation, the Hamiltonian in the "generahzed" 
length gauge without the long wavelength approximation, is given by |46) 



2m 



1 
2m 



ekr 



dAV 



, ^dA dA, ^ 
k-V— + — k- V 

drj arj 



(97) 



The resulting Hamiltonian is similar to the one in the dipole approximation but has three new terms (shown 
on the second line of the last equation). 

The generalization of the Kramers-Henneberger transformation can also be performed. The 
transformation is defined as 

g(all,ia) ._ ( -gQ, . 

where 

1 



A{r]')df]'. 



(98) 
(99) 



This transformation yields a Hamiltonian given by 

1 „2 

2 



A + — A2(,;)+K(r + a) 
Zm Zm 



1 

2m 



da 
k— ■ p 

drj 



dOL 

drj 



p k • p 



(100) 



Again, this gives a Hamiltonian which is similar to the one in the dipole approximation, with some additional 
terms. 

Within this generalization to all orders of the wavelength approximation, the resulting expressions for 
the Hamiltonian is much more intricate and thus, certainly harder to solve. However, they may be useful to 
find other approximation schemes, as discussed in [46j . 



5. Analytical approximations 

There are many approximation schemes that can be used to obtain solutions of the TDSE wave function. 
In this section, we focus on the usual perturbation theory, where the expansion parameter is the weak 
field, or the strong field approximation. We show that observables calculated in different gauges using these 
techniques are invariant under gauge transformations and describe some specific well-known examples. All 
these gauges are simply connected by use of operator equivalent for the total linear momentum in free fields 
[47] p = mv = -ihV = im{Ef - Ei)r/h = ih\IV/{Ef - Ei). 



5.1. Perturbation theory 

In this paragraph, we describe gauge invariance of transition amplitudes when the electromagnetic field 
vanishes at the initial time t = to and final time tf. Although this may seem academic at first sight, this is 
a very important topic because this may be the source of computational errors if proper care is not taken. 
Our discussion starts with the effect of gauge transformations on the energy spectrum and eigenstates of the 
time-independent Schrodinger equation as this is used as a basis for time-dependent perturbation theory. 

Let us consider a general gauge transformation G which relates the wave function in two different gauge 
choices. As usual, the wave functions obey (in this section, we use units in which h = I) 

tdti^p^'^t)) = H^'\t)\^b'^'\t)) , idt\4,^^Ht)) = H^^\t)\4,'^^Ht)). (101) 

Assuming the electromagnetic field vanishes a,t t — t^ and t — tf , we consider that the Hamiltonian will be 
time-independent in these limitt[jj] and thus, the solution of the Schrodinger equation reduces to an eigenvalue 

II This however is not the most general case for the vanishing of the electromagnetic potential. Generally, a gauge where the 
potential is given by A(x, t) = V(/i(x, t) and f7{x, t) = dt<j>(^, t) + C also gives a vanishing potential for ij> an arbitrary function 
and C an arbitrary constant. In this case, the Hamiltonian is not time-independent. 
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problem, for t G {—oo,to] U [i/,oo): 

i^Sl^S)-^Sl<l), (102) 

where the subscript a refers to quantities or operator evaluated at t = to while b is for t = tf (for instance, 

we have H''^\to) = h'^^ and H^^\tf) = H^^). It should be noted here that the eigenenergies for the two 
gauges may differ. This occurs because although the electromagnetic field vanishes in these limits, it is 
possible that the potentials (scalar or vector) still has a non-zero value and this will change the value of the 
eigenenergies: they will be shifted by a certain amount. This can be seen as follows. We have shown earlier 
that the Hamiltonian is not invariant but is changed under a gauge transformation. This transformation 
evaluated at t = tQ,tf is given by 



<i - GHl]lG^ - dtF{^,t)\,^,„^,^ , (104) 

where F is the arbitrary function defining the gauge transformation (G — cxp{iF)). Substituting this 
transformation into Eq. (|103p . we get 

^^Sl-^^S) = [G^Sct - d,F{^,t)\,^,^,^] (105) 

Using the gauge transformation of the wave function and multiplying by l\ on the left, we obtain the 
following important condition: 

AE^,, := Ei'} - Ei'l (^21 (x,t)U,^,^ j^^). (106) 

This expression gives the relation between the eigenenergies expressed in different gauges, in the limit where 
the external electromagnetic field vanishes. At this point, it has been advocated by certain authors that the 
most general transformation should be given by = /(x) (/ here is an arbitrary function of x, independent 
of time) such that dtF = and the energy shift is AE = 0, on the basis that the spectrum should be 
invariant under gauge transformations |48[ I49| . Our point of view on this is different: rather, we assume 
that the term dtF{x, t)\^^^^ is an analytical function of x, such that it can be expanded as a Taylor series. 
We obtain 



oo 



■'^ nx\nJ.nz\ ■ 

where fln^.n^.n^ are the coefficients of the Taylor series and {4'^^b\^'^'"y"^ ^"'I'^^ah 

ria;, riy, n^'th moments. When only the zeroth moment is non zero, every eigenstate is shifted by the same 
amount under a gauge transformation: thus, the latter corresponds to an overall shifting of the zero point 
energy. However, when higher moments are involved, such as the first moment {nx,y,z = 1) which corresponds 
to the average position of the electron in state a, then the energy shift of each state is differen10. Nevertheless, 
the transition amplitudes will be gauge invariant. 

Transition amplitudes between two energy states (i?a^'^\ |0a^'^^)) to (£'^^'^\ l^l^'^^)) are now considered. 

(^121 fl2l 

Physically, this corresponds to preparing the system in states {E\ , |0a )) at i = tg where the external 

fl 21 fl 21 

field vanishes, evolving these states to t = tf and projecting on the states {E^ ,\4>\i ))■ In laser-atoms 

(\ 21 

or laser-molecule interactions, the states cfy^ ^ would typically corresponds to bound states of the atom or 
molecule while the time evolution would include the laser field. Thus, a certain gauge has to be chosen 
to compute these observables. We now show that the transition amplitudes are invariant under gauge 
transformations, which is in fact a particular case of a result that was discussed above. The transition 

^ It may seem counterintuitive that the eigenvalues of the Hamiltonian operator are not gauge invariant. However, we would like 
to stress that these eigenenergies are not observable quantities. Rather, what is observable are the resonances in the transition 
amplitudes: these appear as peak in the spectrum which are detected as spectral lines in spectroscopic measurement. The 
position of these peaks are gauge invariant because they are physical observables. For certain gauge choices, these resonances 
have the same energies as the eigenvalues of H, in which case, it is possible to give a one-to-one correspondence between them 
(as in the length gauge, for instance). This correspondence is possible when the canonical momentum is equal to the mechanical 
momentum, but this is not true in general. 
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amplitude from time to to tf (assuming the electromagnetic field vanishes at t G (—00, to] U [t/,oo)) can be 
written as 

(^(i)(t/)|L^«(t^-,to)|^«(^o)) = {<P^^^\m'Htf,toMi'^), (108) 
A(^) := {^^'\tj)\U^'\tj,toM^'Hto)) = (0p^|C/(^Ht/,io)|0l'^), (109) 

where we take a single eigenstate as initial and final states. Using the gauge transformation on the eigenstates 
and the evolution operator, it is easy to show that A^^^ — A'^'. What is more interesting however is that 
in some cases, it is more convenient to solve the time-independent Schrodinger equation in a specific gauge 
choice, say in gauge 1, but implement the time evolution in gauge 2. Then, we have that 

= {4'^\GHt^'0'Htf,to)G{toM^). (110) 

The gauge operator included in this last expression is very important to preserve the gauge invariance of 
the transition amplitudes in this "mixed" representation and can be omitted only if i^(io,i) = |15[ 124 ) . 
Once we have this expression, it is possible to calculate transition amplitudes using perturbation theory in 
a gauge independent way where the perturbation is the weak external field. However, obtaining a consistent 
perturbation expansion from this starting point requires the resummation of an infinite number of terms, 
which may be impossible in certain cases. This expansion can be performed much more easily if the preceding 
expression is given in the interaction picture, to which we now turn. The latter is required to obtain a 
perturbation theory where each term of the series is at a given order in the expansion parameter. 

It is well-known that the evolution operator in the Schrodinger picture is related to the interaction 
picture by the following relation [50 : 

ui'UtjM e^"'^>^ U^'\tfMe-'"i>, (111) 

t>l'i,.(i/,io) = e^^">^^(^Ht/,to)e-^^-*«. (112) 

where the interaction Hamiltonian appearing in the evolution operator ([//(O, t) = Texp{J^ dTHi{T)) with T 

being the time-ordered operator) are given by Hj^^l{t) = H^^''^\t) — H^lj^\ When we have Hji''^^ = H^^'^\ 
the transition amplitudes in the interaction picture simply become 

^(1) e-(^^^'*/-^i^'*o)(^«|C/j^)(t/,to)|0i^)), (113) 

^(2) ^-.K^>.,-^(^>.o)^^(2)|j^(2)^^^, ^^^|^(,)^_ (^^4^ 

In the more general case, where ^ H'^'', first, we need to use the property of the evolution operator 

that U'^^'^\tf,to) = U^^''^\tf,t)U^^''^\t,to) for an arbitrary time t &]to,tf[. Then, using the transformation 
to the interaction picture, we get 

A(i):=e-(^^^'*/-^i^'*«)(0«|C/«(t,,t)e-^^^'*e^^^^^ (115) 

^(2) e~^^^^'^'^^-^'^'^^'\<^f\uf}^^^ (116) 

Now, as in the Schrodinger picture, we write the transition amplitude in the second gauge choice as 

A^^):=e-^^^'^''^-^'^'''^^\<^'i^\G\tf)U^^^^^^ (117) 

This is the most general formula allowing to derive a perturbation expansion which is obtained by expanding 
the evolution operators and the exponentials. It is very important to note that if the gauge operator contains 
the electromagnetic potential, G should also be expanded to obtain the "true" leading order contribution. 



5.1.1. Transition amplitude: a specific example Although the transition amplitudes are identical in theory, 
some important differences may occur between different gauges when they are evaluated explicitly. This 
occurs because this usually requires the evaluation of infinite sums on intermediate states and although 
these sums yields the same result, their convergence depends on the gauge chosen. We will show this feature 
in a specific example which involves the the one and two-photon processes. We compare the perturbative 
calculation in the length and velocity gauges. As they are related by a unitary transformation, their 
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corresponding transition matrices are then identical at resonance or not. Again, we refer to [TS^ for more 
complete explanations and calculations. 

Following |15) . we consider a linearly polarized laser pulse in the dipole approximation for which the 
vector potential, in the velocity gauge, can be written as A{t) = A{t) cos{ujt)ex where t e [0, T], with Tuj » 1. 
Here, A(t) is an envelope function which gives the temporal pulse shape. We also assume the presence of 
an atom or molecule having bound states and characterized by a time-independent scalar potential (pc- We 
recall that the transition element from {(pa) to is written as [15] 

5a.6 = ^ii'"^"' = lim hm (0i''"^'^'|C/^''''^°^(^2,ti)|0^"'"^), (118) 

where is the evolution operator for (l) length, (v) velocity and (a) acceleration gauges. It is very 

convenient to evaluate the eigenstates of the system in the length gauge: in this case, the vector and scalar 
potential at t = ±cx) vanish for any pulse shapes because the electric field in these limits is zero and thus, 
the eigenenergies correspond to the physical bound states energies. Then, according to the discussion of the 
last section, we can write [HI [24| 



Sa,b= lim lim {4'^\U^'Ht2,hM^), (119) 

= lim lim (0(')|G("')t(i2)[/W(t2,ti)G(''')(ti)|0W), (120) 

= lim lim {cpi''>\G<^'^'^Ht2)U^''Ht2,h)G<^'-HhM^). (121) 

Then, these expressions can be expanded in powers of e\A\ to get an approximation of the transition 
amplitudes. Here however, it may be more convenient to use the interaction picture. It should be stressed 
again that since G'^^"^ and G'^'"-' contain the potential, they should also be expanded. 

It can be proven (see jl5j), from a perturbative approach at order 1 in e (for 1-photon process) and 
order 2 in e (for 2-photon processes) and under the dipole approximation, that (here and in the following, 
we define hujab = Ea — Ei,): 

• for a resonant 1-photon process sj^^^ = sj^^ ~ ^!ib ■ 

• for a resonant 2-photon process, the transition element satisfy: 

.H_.22^n(.)(^)^,(,^^_2,), (122) 



where the infinite sum over intermediate states is given by 

_ ^ {M^x ■ Y>\(l)r){<Pr\ex ■ P|0a) , . 

and where are the transition states. 

In the case of the length gauge we obtain an equation of the same form: 



^S=e^|QS(^)'^(-.a-2.), (124) 



but now, the sum is 

= 4^ • (12^) 

Using the fact that for any \(t>s)^ \4>t), we have 

{(psWx ■ v\<l^t) = ii^stm{(t>s\ex -Acpt) , (126) 

we can deduce 

^C.) ^ ^2 ^ Jbr^ar ^^^^^^ v\<t>r) {<t>r\Bx ' v\<t>a) ■ (127) 
— n[UJ — LOra) 
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Out of resonance, when 2lj ^ uJab, the transition amphtudes are equal 5*^2 = ^lib ~ because the delta 
function 5{lu — ujab) has a support only at resonance. It can then be proven that the contributions 
and Q^^^ are equal as infinite series at resonance, when 2a; = uiab |15| . Therefore, clearly, the transition 
amplitudes are gauge invariant to second order in perturbation theory in the example considered. 
However, to obtain an explicit result for an observables related to these transition amplitudes, only 
a limited number of intermediate state can be evaluated: the sums has to be truncated because they 
cannot be evaluated in closed-form. It is clear that the convergence rate of series Qib and Qlh'' are 
dependent on w. That is: for a; <C 1, in the length gauge, we notice that the series converges more 
rapidly than in the velocity gauge as is (much) smaller than Wfcr'^ra, for energy states \(j)r) such that 
— £^a| ^ 1 and \Er — Eb\ ^ 1. The opposite occurs for a; ^ 1. Denoting 

^abr = A (t>a\ex ' r | (/)^ ) ( (/)^ | • r|(/)a), (128) 

h[u - UJra) 

Up- 

Qal = TT -Me. ■ r\(br){A\ex ■ r|</),), (129) 

n[U! — UJra) 

we get the following results: 

— when cj is a low frequency compared to the transition energies ujab, for most r 

4t= «!' (130) 

g^"") UJbr^^ar 
^abr 

SO that length gauge has a faster convergence. 

— when w is a high frequency larger than ujab, for most r 

'iabr 

SO that velocity gauge has a faster convergence. 
This justifies the use of distinct gauges for different regimes. 

5.1.2. Dynamic Stark shifts In the following we shortly focus on the frequency dependent shift in energy of 
any nth atomic level, see [51] for details. First we recall that the elements of the radiative matrix, for states 
\4>a), \<t>b) are written as: 

ih=— = {cj)a\[v,HQ]\cj)b) = hhJab^ab, (132) 

m 

where hjjab = Ea — Eb and Hq the laser-free Hamiltonian. This gives directly a relation between the 
momentum radiative matrix elements (velocity gauge) and the dipole radiative matrix elements (length 
gauge) and A{t) = Aq cos(wt)e: 

• In the length gauge: EoVab- 

• In the velocity gauge: 

=z r^b-E, (133) 

mc UJ 

The equality with the length gauge, occurs only at resonance uj — uJba- This gauge is then appropriate 

for Uj/uJba < 1. 

• In the acceleration gauge for perturbative fields 

V{a{t) + r) = V{v) + a{t) ■ VV + (134) 
with ih — W = [i?o,p], then (Vy)afc = i^baY^ab, we obtain 

a{t) ■ {^V)ab = l — i—Vab) = -"^E^Vab. (135) 



mcV UJ / uj'^ 

Again, looking at near resonance, there is equality with the 2 other gauges. We then deduce that this 
gauge provides an even faster convergence than the velocity one for uj/uJba *C 1 
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The corresponding energy shifts, in the length gauge are (see [5T]): 

2 ^ LOnm ± 



(Ai?l'^)« = ^E^^|Eo-r™„P, (136) 



where hunm — En — Em- For large laser frequency uj, the energy shift tends to the ponderomotive energy 
Up which is positive. That is for uj ^ ujnm 

(A£;^'))„ ~ -_^,^„„|Eo •r„„|2 = AE, ^= C/p. (137) 

m 

This is a simple consequence of the Thomas-Reiche-Kuhn sum formula ~ ~l/2- The total 

energy shift contains a continuum, /, contribution and a bound state, b, contribution 

AEc = A£;^'^ + A^y^ . (138) 
from which we deduce that 

{AEf\ = -^Y: ^J5%;i|Eo • (139) 

and for large frequencies 3> 1 

Similar computations in the acceleration gauge, leads to, for a; ^ 1, 

(^^i^^^ - ■ r„™P - 4.gp„(0), (141) 

where pn denotes the nucleus electronic density and explain the Lamb shift of s atomic states. 



5.2. Strong Field Approximation (SFA) 

SFA is a very common non-perturbative approach to study the interaction of bound systems with intense 
lasers: only the ground state and continuum contributions are taken into account to represent the 
wave function for an intense laser-molecule interaction and is usually based on a Single Active Electron 
approximation, which can be extended to multi-electron systems Early work on the equivalent gauges 
for strong field physics has emphasized the difference between the length and the velocity gauges. It was 
shown that the TDSE expressed in center of mass (cm.) coordinates and relative coordinates is nonseparable 
beyond the dipole approximation in the length gauge, whereas in the velocity gauge, separation of variables is 
straightforward [52 [ l53 j . In the dipole approximation and length gauge, the radiative interaction is described 
by a scalar potential — E • r whereas in the velocity gauge the interaction involves a gradient, A • V, 
reminiscent of nonradiative interaction in molecular physics [2]. Recently the acceleration representation 
has been generalized beyond the dipole approximation 1541 as this representation is most convenient for high 
frequencies where "atomic stabilization" is expected to occur at high intensities. The SFA theory of HHG 
is basically a time-dependent two potential Distorted Wave Born- Approximation (DWBA), where the final 
state of the electron can be described as a time dependent dressed state, called a Volkov state 55 , whereas 
the initial state, due to its high ionization potential energy /p, is considered unperturbed. Calculations of 
photoionization based on Volkov final state wave functions were first done by Keldysh [56], and later by 
Faisal [57l and Reiss [58^ . It was also re-examined recently for circular polarization [59l |60] . This strong 
field approach produces very accurate results for multiphoton detachment for negative ions since Coulomb 
potentials are absent in the ionization process in both length and velocity gauges [6T1[62] . Later work by Bauer 
et al. |63| showed that the prediction of the two gauges can differ qualitatively for ionization of negative ions 
and that the length gauge SFA matches the exact TDSE numerical solutions J33| . More recently, it has been 
shown that the SFA used in models for computing HHG invalidates the Ehrenfest theorem [64] . Gordon and 
Kartner have shown in fact that the SFA can be improved by using the acceleration radiative interaction W 
in the emission matrix element since then the SFA HHG amplitude is correct to first order in the Coulomb 
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potential V{r) [SS]. The SFA encounters further difficulties in interpreting Molecular High Order Harmonic 
Generation (MHOHG) spectra due to multi-center nuclear interference effects both in the ionization and 
the recombination processes |66[ 167] . In the acceleration representation, each nucleus contributes to the 
MHOHG amplitude via the multicenter nature of the total Coulomb force, —WV. This multicenter effect is 
absent in both length and velocity gauges [6 7) . 

The non-invariance of SFA under gauge transformation is an issue which is also studied in [6 8) , |69| . We 
here give some details about this issue. As often noticed (see also [70]) and recalled above, the ionization rates 
and energy distributions of strong-field photo-ionization differ depending of the choice of the gauge (length, 
velocity or acceleration). However, there is no definitive conclusion regarding the most appropriate (all give 
relevant results depending of the physical problem under consideration [7T], [53]). We note that Faisal [72] 
showed that length and velocity gauge SFA amplitudes under the dipole approximation are equivalent at 
all orders under some assumptions on the initial and final state partitions of the Hamiltonian. Details of 
what follows can be found in [71] ■ We recall that the ionization amplitude from |'i/'o) (ground state) to the 
continuum jf/'/) is defined as 

So,i= lim lim {<j)i\U{t2,h)\(bo) . (142) 
The propagator [/ is a priori gauge dependent and is defined as follows: 

U{h,t2) = Uoih,t2)~-i f dTU{ti,T)Hiaser{T)Uo{T,t2), (143) 

where Uq is the time-evolution propagator associated to the laser-free operator and H^^^^^ is the interaction 
Hamiltonian that includes the laser field in the length or velocity gauge, respectively. Then, the transition 
amplitude is given by 

^0./ = -* lim r dT{cj,j{t2)\U{t2,T)Htl{r)\Mr))- (144) 

The SFA approximation for the ionization amplitude consists of replacing the exact final state wave function 
\ipi} by the Volkov one {ipvoikov)- 

{lJjl\U{ti,t2) ^ti^-oo {lpvolkov{t2)\, (145) 

with A(±oo) ~ 0. Under some assumptions [71) 's authors show that the two representations (in fact in any 
gauge) are equivalent to the length gauge representation. 

So.i = -i I dt{4,,oikov{T)\H'i'^^l{T)m)). (146) 

As mentioned in [63 , for t < t the electron is bound and the laser-electron interaction is neglected. For 
t = T, ionization occurs and the electron then moves rapidly out of the range of the Coulomb potential. 
Keitel et al. have studied a gauge-invariant relativistic version of SFA in [75]. They also provide a gauge 
invariant ionization amplitude expression which coincides with the SFA in the length gauge, considered as 
more adequate for ionization process (ATI or nonsequential double-ionization (NSDI)). 

6. Numerical approximations 

This section is devoted to the numerical computation of TDSE in the velocity, length and acceleration 
gauges. The motivation comes from the fact that although in theory, the transition probability, dipole 
moments, velocities and accelerations are equal in all the discussed gauges, in practice due to approximations 
(continuous or discrete) the equality of these quantities does not hold in general. In addition, the choice of the 
mathematical structure (which is gauge dependent) of the TDSE guides the choice of the numerical method. 
TDSE in the velocity gauge, contains a transport operator which is more accurately approximated using finite 
difference or finite volume techniques. In contrast, the Laplace operator is particularly well approximated 
using variational techniques, such as finite element or spectral methods. These aspects will be discussed in 
this section. For each approach, the principle of the numerical method is summarized and an appropriate 
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application (at least the first steps of the discretization) of TDSE's is proposed. Before detailing numerical 
methods however, we describe the non-perturbative process of harmonic generation. This phenomenon, 
which is usually studied theoretically by using numerical solutions of TDSE, is an example where the gauge 
choice is instrumental in obtaining accurate results. 

6.1. Harmonic generation 

In the following, we discuss the dipole moment, velocity and acceleration calculation in the dipole 
approximation for HHG processes. Prior to this we recall the Ehrenfest's theorem. For any operator O 

\{0)^^{[d,H]) + l"—). (147) 




dt ih 

In the dipole approximation, one can derive simply the electric field E{t) via the vector potential: 

A{t) = -c£{t) sin {u{t-tc))/uj. (148) 

and 

E{t) = ~dtA{t)/c = £{t) cos {uj{t ~ tc) +4>) + Ecor , (149) 

where tc is the peak of the field and cf) is the Carrier Envelope Phase (CEP), £{t) is the envelope of the pulse 
and 

Ecor = sin {uj{t - tc) + <j))dt£{t)/uj, (150) 

is the correction to the simple form (I149p arising from the derivative of the envelope of the vector potential. 
Ecor is negligible near the peak, t — t^ for long pulses and for slowly varying envelopes £{t). For short pulses 
such that the envelope varies rapidly, one must use the complete form (jl49p . which ensures the zero-area 
theorem: 

E{t)dt ^ A{ti) ~ A{t2) ^ . (151) 

tl 

One usually chooses A{ti) = A{t2) = that is £{ti) — £{t2) =0 in (|148|) . If the total area of the electric field 
is not zero, a simple Fourier transform of (jl5ip shows that Jj^ E{t)dt — Eq, that is a static field component 
is present in the pulse, contrary to Maxwell's equations. 
In general the dipole velocity is related to the momentum as 

z = -i[z,H]=dH/dp^^p^. (152) 

for which we define {z{t)) ~ ('i/;(t)|p^|?/'(t)) andp^ ~ —id/dz. The Fourier transform i 
is given by 

z{uj) ^ e''^\z{tf)) +iLoz{Lo) (153) 

where the initial condition {z{ti)) = for symmetric states. A similar procedure for the acceleration |74] 
gives 

-£ = -i[z,H] = -i[p^,H]= ~dH/dz, (154) 

then 

z{io) = zcji(cj) - e'"*^ (i(i/)), (155) 

with the initial symmetry dictated condition {z{ti)) = 0. We note that this condition differs from that 
imposed in the tunnel ionization model |101 lll[ I74[ I75| . which assumes instantaneous initial condition 
z{ti) = 0. Equations (|150p . (|152p demonstrate that the three forms of the particle operators involved in 
radiative interactions, z{uj), z{uj) and z{uj) may differ radically depending on the final values of the average 
dipole {z{tf)) and velocity {z{tf)). This issue originally raised by Burnett [75] was discussed in detail in [TJ] 
and we show another example in Figs 16.11 16.11 how these differ significantly on the parameters that specify 
the laser pulse (duration, intensity, CEP). Figures [01 and [01 illustrate the effect of pulse duration on the 
X and y components of MHOHG from the one electron molecular ion at internuclear distance i? = 22 
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a.u. ionized by a circularly polarized pulse at wavelenth 400 nm(?ii — 0.114 a.u.), and intensity I — 10^'^ 
W-cm~^ incident on the x, y plane of the molecule with the x direction parallel to the intermolular axis, i.e. 
R. The internuclear distance R = 22 a.u. is chosen to coincide with the laser induced collision radius for an 
electron ionized from one end of the molecule with its neighbor as predicted by a classical collision model. 
Px^y{uj) correspond to the FT of the squared average value of the dipole, velocity and acceleration operators. 
We note in Fig. 16.11 for an ultrashort 2 cycle pulse that the MHOHG spectrum is the same for dipole and 
velocity since both are nonvanishing after the pulse whereas the acceleration which does converge to zero 
average value gives a completely different spectrum. Fig. 16.11 corresponds to a long 20 cycle pulse for which 
all three average values of operators vanish after the pulse, thus giving basicallly identical MHOHG spectra. 
Figures ISTTl 16.11 confirm the numerical difference in HHG calculations in different gauges for ultrashort pulses 
due to the nonvanishing of dipoles or currents after the pulse. 

If the final zero condition, {z{tf)) = {z{tf)) = are satisfied, then the three forms of HHG spectrum 
intensities corresponding to dipole, velocity and acceleration are simply related by 

Di{uj) ^ uj^Diiu;) = uj^D^iuj), (156) 

where D^{uj) — and £_{uj) ~ J^^ (iie^'"*^(i). Earlier work on calculations of tunnelling rates by 

TDSE's or time-independent methods emphasized the difference in switching E(t) or A{t) and obtaining 
gauge invariance tunnelling rates |77J [7S] . 

A more recent discussion of the related problem of adiabatic cut-offs of fields in gauge invariant 
electrodynamics emphasizes the importance of the electromagnetic potentials (|148l) , to generate the physical 
field (UMl) l^- 

Finally as pointed out in the previous sections, the dipole and velocity forms of radiative interactions are 
considered as gauges as they give rise to gauge invariant Lagrangians, whereas the acceleration form is 
generated by a unitary transformation in TDSE's, and is therefore not a gauge transformation. Recent 
work has emphasized the direct relation of the velocity rather than dipole or acceleration to the harmonic 
fields generated in HHG [HOI ISB IH2]- Numerical simulations comparing the scaling in HHG intensities 
in the velocity versus the dipole gauge have concluded a more favorable scaling in the velocity gauge by 
noting that the canonical momentum is reduced by the vector potential since v — p — eA/mc [30 . On 
the other hand, as already noted by Gordon et al. [SS], the acceleration representation of HHG provides 
considerable computational advantages in SFA as the representation includes Coulomb potentials to first 
order, which is absent in the two other gauges. Similarly in MHOHG, multicenter Coulomb effects are 
essential to predict maximum in the spectra and these are correctly predicted by the acceleration form 
|39i 166] ■ Furthermore, the acceleration representation offers the computational advantage in non Born- 
Oppenheimer molecular simulations, that contrary to the divergent length gauge radiative term zE{t), 
radiative terms vanish asymptotically so that projections can be made simply onto free electron Coulomb 
states |83j . Comparison of convergence between the same discretization schemes for all three representations, 
dipole, velocity, and acceleration have demonstrated the superiority of the acceleration form due to its 
similarity to Lagrangian adaptive grid methods used in fluid dynamics |84j |85] . 

In |30| . the authors compare numerically the dipole moment computed using the length and velocity 
gauges for the hydrogen atom. As expected a very good agreement is obtained in these 2 gauges: 

i^^ir,t) = -(A + -—^-——)^ir,t), (157) 



dr^ ' ' \ \v\+W{Y,t) 
with in the length gauge (linear polarized laser field) 

VF(r,i) =r-E(t), (158) 
and in the velocity gauge 

PF(r,i) = -iA(t) • V(t), (159) 
where the laser pulse is chosen as: 

E(i) - i;o(sin2 (^) sin(L^ot) - ^sin (^cos(L^oO))e., (160) 
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Figure 1. Dipole, velocity, acceleration for a 2 cycle circularly polarized pulse at 400 nm and / = 
10^''W-cm~^ and corresponding HHG spectrum for H2 at distance R = 22 a.u. 



with E,{t)dt = and 

A(t) = ^sin^ {%) cos(c^oO))e., (161) 

A = 800nm and T — 110. 32a. u. and / = 3x lO^^W-cm^^. The equation was rewritten in spherical coordinates 
and the chosen numerical method is was spectral/finite clement method with B-splines basis for the radial 
part and spherical harmonics for the angular part. Due to the choice of linear polarized laser field, the 
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Figure 2. Dipolo, velocity, acceleration for a 20 cycle circularly polarized pulse at 400 nm and I = 
10^*W-cm~^ and corresponding HHG spectrum for //^ at distance R = 22 a.u. 

azimutal magnetic quantum number was taken null. Let us shortly recall the conclusion of this interesting 
paper. To obtain convergence (up to a certain fixed error): 

• More grid points were required in the velocity gauge and larger angular basis in the length gauge. 

• At high intensity the convergence was faster in the velocity gauge than in the length gauge. 
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Notice that the autors expect for nonhnearly polarized laser fields the computations in the velocity gauges to 
become more striking. Naturally, the convergence as function of physical parameters, is also highly dependent 
on the choice of the numerical method to solve the TDSE. The above conclusions are of course a priori not 
valid for other methods, see [55] . 

Note that in numerical computations also show that the acceleration and dipole forms of the 
transition matrix lead to different results. It is shown that the acceleration form accurately predicts the 
interference effects, exp (ik • i?/2) ± exp ( — ik • -R/2) in iJ^, while the dipole form does not (see table II 
and eqs (34, 43) in [33] )■ This is due to the fact that the acceleration operator is made up of the forces 
from each Coulomb center, thus representing to first order the multicenter molecular structure, whereas the 
dipole operator emphasizes larger distances and contains no molecular structure information. In addition, 
the dipole form of the transition matrix indicates that the emission of harmonic photons is proportional to 
"momentum" derivative of the Fourier transform of the ground-state wave function. This feature is a key 
ingredient of the method for tomographic imaging of the ground-state molecular orbitals presented in [9] by 
exploiting the orientation dependence of MHOHG spectra. According to [33] a formulation of tomographic 
imaging based on the acceleration could then be more accurate in reconstructing molecule orbitals from the 
harmonic spectrum. 

When dealing with numerical computations of transition amplitudes, it is important to consider 
sufficiently large computational times. Indeed, theoretical estimation of these amplitudes necessitate in 
theory a computational time from — oo to -foo. Time truncation introduces errors. For this reason in [87], 
[55] the authors focus on the solving of TDSE's using the acceleration gauge. 



6.2. Numerical methods 

We start this section with a discussion of a fundamental tool used to numerically solve PDE, the time 
splitting method. We then present several numerical techniques which are commonly used to solve TDSE's 
in different gauges. 



6.2.1. Splitting The TDSE involves the coupling of several operators of different type: hyperbolic 
(9t — A • V), parabolic {dt — A) and algebraic {Vc or r • E). Direct (no splitting) approaches are commonly 
used and will be discussed later, however solving individually these equations can also lead to very accurate 
numerical solutions (spectral convergence, exact). We start by detailing the general principle of splitting 
methods. 

• Mathematical Splitting: Trotter's formula [89]. We consider the following equation: 

M 

ut^^AkU,u{r,0)^uo{r) (162) 
fc=i 

where {Ak)k={i^--- ,m} is a finite sequence of (spatial differential and algebraic) operators. Formally 
the solution (via time propagators) can be written as u{r,t) = e* ^'==1 '^'=uo(r) if the operators are 
time-independenQ The exponential is defined by the Trotter-Kato's formula which states that 

e*E£,A. ^ fn^,e*^'=/^) . (163) 

• Numerical Splitting. A discrete version of Trotter-Kato's formula is the well-known time-splitting 
method. From [0,dt] and rather than solving (|162p [55] . we set 

(164) 



Ut = 


AlU, 


M(r,0) 




Ut = 


A2U, 


w(r,0) 


= ui{r), 


Ut = 


Amu, 


' ' ' 1 

M(r,0) 


= UM{r) 



+ If the operators are time-dependent and do not commute at different times, the solution will be given by ii(r,t) = 
Texp I^/q cii' X]ft=i ^fc "o(i")i where T is the time-ordering operator. Splitting methods can be easily extended to this 
latter case. 
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where Uk{r,dt) is solution on [0,(it], to 

ut^AkU, u(r,0) =Ufc_i(r), fc e {1, • • • ,Af}. (165) 
This in fact consists of applying Trotter-Kato's formula as follows: 

edtj:tUA,^ lij^^ (u^I^c'^'^./nY , (166) 



and taking = 1. The error between the exact to (I162p and approximate to (|164l) solutions can easily 
be evaluated and is equal in that case to: 

u{r,dt) = e'**^"i^*wo(r) = Ilf:l,e'^*^^uo{r) + 0{dt^). (167) 
Errors are due to the non-commutativity of operators, [85 . More precisely for M = 2 

i Ut^Aiu, u(r,0) = uo(r), , 
I Ut^Mu, u(r, 0) = ui(r), 

For t e [Ojdt], we formally have the identity 

exp {dt{Ai + A2)) = exp(dMi) exp(diA2) + A2]/2 H , (169) 

where [Ai, A2] — A1A2 — A2A1 is the commutator. When Ai and A2 commute the splitting is exact. 
More generally for M operators, if [Ak, Ai] = for all k and I in {1, • • • , M}, splitting the equation in 
M equations does not introduce any error. 

Higher order splitting approaches are naturally possible, such as the famous Strang splitting [90] , 
Considering the equation 

Ut{x,t) = Aiu{x,t) + A2u{x,t) on [0,dt], u(r, 0)=uo, (170) 

where Ai and A2 are two algrebraic or differential spatial operators, the principle of Strang's splitting 
consists of approximating the exact solution Uoxact(r, dt) = exp ((Ai + A2)dt^uo{v), by solving 

ut = Aiu on [0,dt/2], u{r,0)—uo{r), 

ut = A2U on[0,dt], u(r,0) = ui(r), (171) 
Ut — Aiu on [0,dt/2], u{r,0) — U2{r), 

The calculated solution is given by 

-"approximate (r, dt) = cxp {Aidt/2) cxp (^2^^) cxp [Aidt /2)un{r) . (172) 

It can easily be proven that in that case the error between the exact and approximation solution is a 
0{dt^). In practice, we will take Ai = A, A2 = Vc + A • V or A2 = K + r • E. The advantage of this 
approach comes from the fact that ut = A^u (fc — 1,2) can be very accurately solved, at least more 
than Ut = Aiu + A2U. The price to pay is that using a splitting method requires to choose dt sufficiently 
small to reduce the error. The splitting error is obviously added to the discretization errors coming from 
the approximation of equations (|164p . Splittings are extensively used with spectral, real space or even 
exact methods. 

The generalization to high order splitting methods for solving TDSE's is described in [91], [38] . 

Before describing other numerical methods and discretization, we propose a short discussion which attempts 
to link the TDSE to fluid dynamics equations (such as the transport and Navier-Stokes equations). These 
fluid equations have been studied extensively by mathematicians and physicians, and many techniques have 
been developed to analyze and solve numerically these equations. Thus, the formal analogy existing between 
TDSE's and fluid dynamics equations allows to use these techniques in the context of quantum mechanics, 
[84l |45] . Our remark is more specifically devoted to find connections between the TDSE written in the 
velocity gauge, with usual transport problems as well as incompressible viscous or non-viscous fluid flow 
equations. The TDSE involves a kinetic operator. A, as well as a transport operator V: 

i—^p^-—Ai: + iA-Vi: + —A^'th. (173) 
ot 2m 2m 
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It is possible to solve numerically this equation splitting it into two parts (following the previous discussion): 
4^.-i-A^ + i-AV, (174) 

and 

^ + A-VV' = 0. (175) 

Under the Coulomb gauge, the last equation is a transport equation equivalent to the following conservation 
law: 

d 

+ div(AV') = 0. (176) 

The general solution to this equation is obtained using the method of characteristics. Defining the 
characteristics as integral curves of 

X(t)=A(X(i),i),X(0) = Xo, (177) 

along these curves the solution is constant — (X(i),tj — 0, which allows to deduce the exact solution. This 

approach naturally fails when viscosity (real or complex) is added. 

Interesting connections can be found with incompressible fluid flows. Under incompressibility condition 
V • u 0, the conservation of momentum of a non- viscous fluid of density p, becomes 

d 

— u + div(pu ®u + P) = Q, (178) 

where P is the pressure and u the fluid velocity. Although it is tempting to generalize the comparison 
between TDSE in the velocity gauge and the momentum equation for incompressible flows, the solution 
type for p76p (which is linear) and (|178p (which is nonlinear, in fact quasi- linear) maybe very different |92) . 
Including now viscous effects, the Navier-Stokes equations is written as: 

p—vL + pu • Vu = V • cr, (179) 

where for Newtonian fluids, the fluid viscosity /i is constant and V • a writes /xAu — VP. For non-newtonian 
fluids /i is no more constant and can even be complex. We refer to ;93. for the interested readers. It is then 
interesting to notice the presence of a "complex viscous" term, which allows to make a direct connection 
via the complex viscous and transport terms, with ()173|) and ()174|) . Therefore, some mathematical and 
numerical techniques to solve the TDSE can be derived or adapted from fluid mechanics. The interested 
readers could explore further this question. 



6.2.2. Galerkin Methods Spectral Methods. The principle is to search for the wave function ip in the 
form: 

./-(r,*) =^c„(O0„(r), (180) 

n 

where ((/>„)« is a basis of L'^{R^) containing the exact wave function. Note that by doing this, the continuous 
states become discrete in the numerical calculation. The Galerkin methods are more adapted to solve TDSE 
in the length gauge than in the velocity gauge. This comes from the fact that for stability reasons, the 
discretization of transport operators (such as A • V in TDSE written in the velocity gauge) necessitates a 
numerical upwindind (see details in Section r6.2.3|) . By default, Galerkin's methods are centred and as a 
consequence provide unstable approximate transport operators (some complex empirical techniques exist, 
such as the Streamline Upwind Petrov Galerkin method). The Galerkin method proceeds in the following 
way. We formally write the TDSE as (i? = - A + K + r • E) 

= i?V'. (181) 
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We multiply by a test function (j>k G i^(IR^) and integrate on space, which transforms the partial differential 
equation into a variational formulation. Denoting (•, •) the inner product in L^(R'^), that is 



0^(x)0i(x)d3x, / |0|2(x)d3x= (182) 



this becomes 



iVc„(O(0n,</)fc) = > 7770„,0fc), (183) 



which becomes, by integration by parts, and assuming that the basis functions vanish at infinity: 

*EnCn(i)(0n,</>fe) = - E„ C„ (i) ( V(/)„ , V(?!)fc) 



(184) 



+ E„c„(t)([K(r)+r.E(t)]0„,0fc). 

Trucating the sum and keeping the N first terms, the set of equations can be rewritten as: 

AC{t) ^{B + D{t))C{t), (185) 

where A = ((0„, 0fc))„^fc, S = ((V0„, V(^fc)(i2(R3))«)„ and D{t) = ((14(r) + r • E(t)</>„, </)fc))„ are Afjv(C) 

matrices and C{t) = (ci,{t),--- ,CN{t))'^ G C^. The time discretization leads to the solving of a linear 
system. For instance if a simple forward Euler discretization of the time derivative is applied, (|185p becomes 

AC"+^ ^ AC" + dt{B + D'')a\ (186) 

More elaborated time discretizations are of course more appropriate from a stability as well as accuracy 
points of view. Sparsity of the matrix A (many zero entries) is crucial for computational efficiency (data 
storage and computational complexity). What characterizes spectral methods is the non- locality of the basis 
function, which have support in all the spatial domain. Different choices of basis are possible, the most 
common are: 

• '/'n(r) — exp(inr), corresponding to Fourier series expansion. 

• all kinds of orthogonal polynomials (Hermite, Legendre,...). 

• Spherical harmonics. 

• Eigenfunction decomposition. This method requires the eigenfunctions: 

H0(t>n = Er,(t>n. (187) 

Although this is the most accurate basis (as the solution is decomposed on the exact orbitals), it 
requires important preliminary work, consisting of determining large sets of bound states of the laser- 
free Hamiltonian. Perturbation theory for solving TDSE's is also based on this decomposition. 

The spectral approach is particularly appropriate for approximating the kinetic operator, in particular 
using a Fourier decomposition (which transforms the kinetic operator into -|jk|p). Spectral (or exponential) 
convergence is possible in general, and Gibbs' phenomena (oscillations near singularities) do not appear (in 
general) due to the regularity of the solution. The consequence is then a fast convergence. For any smooth 
function ip, say 27r-periodic, a A'^— term Fourier series approximation is such that (spectral convergence) 

II V - II C{q) exp(-A)|| V|U2([o,2.]), (188) 



where 



and 



i'hi^)^ ?A„exp(mx), (189) 



2-77 



ipn = — / i^ix) e'x:p(—inx)dx. (190) 
27r Jo 

These methods are then very attractive (easy to implement and very fast convergence). We refer to [94] for 
interested readers. 

31 



Finite Element Methods (FEM) [95], [96], [97]. As a Galerkin method, the finite clement approxi- 
mation is very similar to spectral methods. The main diff'erence comes from the fact that the basis functions 
have a bounded support and are only piecewise regular. Again, this method allows to consider non-regular 
domains and non-uniform meshes (useful to capture singularities or large gradients). 

*E„c„(t)(0„,(?!)fc) = -E„c„(i)(V(/«„, V(?!)fc) 

(191) 

+ E„c„(t)(T4(r)-Hr.E(i)<^„,<^fc) 

the basis functions {(j)k)k are for instance, piecewise Lagrange polynomials (or B-splines [551l5^[l00l ) with a 
localized support. This allows in particular to increase the sparsity of the "mass" ((/ 4>i(t>j)ij) and "stiffness" 
{{V(j)i\/(j)j)ij) matrices. In addition non regular solutions are more accurately captured compared to spectral 
methods. Many convergence results exist, in particular for Lagrange finite element methods. The basis 
functions are piecewise polynomials of degree k equal to 1, at one node and otherwise. Typical error 
estimates for the wave function in Sobolev space, that is for ^ e H''{R'^) (V«V e {L^iR"")) , I = 0,...,fc) 
and piecewise approximation by polynomials of degree k is: 

U-M^Ch''+'UU2^[o.2.]), (192) 
for h, largest cell area and C a positive constant. 

Collocation Methods. Starting again from 

i^(r,t) = ^c„(i)0„(r), (193) 

n 

where ((/>n)n is a basis of L^(M'^). We again formally write the equation 

= H^p. (194) 

We multiply this time by test functions which are (5-functions in freely (typically finer near singularities) 
selected grid points (r^): that is defining by 6k ~ S{rk), we get 

*E„c„(i)(0„,^?fc) = -E„c„(O(A0„,0fc)-f E„c„(t)(K(r) 

(195) 

+r •E(i)0„,0fc) 
which this time leads to a set of equations: 

i^c„(t)0„(rfc) = -^c„(OA0„(rfe) -F^c„(0(Fc(rfe) +rfc •E(0)<^„(rfc). (196) 

n n n 

In particular a space discretization (on a unstructured grid) of the Laplacian is then necessary. This is 
usually done using Taylor expansion techniques. The time discretization leads to a system of the form (jl86p . 

In general, these approaches do not necessitate any splitting to be highly accurate. Due to the total 
freedom on the grid point locations, they are very attractive for approximating (smoothly) singularities. The 
full convergence analysis of this approach is however still largely open. 

6.2.3. Direct Real Space Methods Finite Difference Methods (EDM). These methods constitute the 
most simple approaches to discretize TDSE's. The principle consists of approximating the spatial derivatives 
as follows. In Tq — (xq, yo, zq) and time to 

a ^ , ip{xo,yo,zo,to) -ip{xo- dx,y(),za,ta) 

OxV[Xo,yo,zo,to) ^ J , (197) 

ax 

or 

„ ^ , ilj{xo + dx,yo,zo,to)-tp{xo-dx,yo,zo,to) 

OxV[Xo, 2/0, zo, to) ~ T-j . (198) 

2ax 
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Consistency (correct approximation of the equation), stability (numerical solution remains bounded) and 
accuracy questions are discussed in details in [lOlj . Typically, Crank- Nicolson's scheme is used to 
approximate (here in 1-d) TDSE in the length gauge: 

dti^ = td^^^-iVc{x)^P + ixE{t)^P, (199) 

Denoting by "0" an approximation of •ip{xj,tn) {xj — jdx, t„ = ndt), the scheme writes 



^(^r^-^.") = -i(x,i?(Wi) + K(x))0r' " 2lx^(^7ti-2i^^'+^^i), 



(200) 



This scheme is of order 2 (error divided by 4 when space step is divided by 2) in space and time and is 
unconditionally stable (that is, is stable for any choice of da; and dt). Stability is an important criterium with 
consistency (that is the numerical scheme approximates the continuous equation) to ensure the convergence 
of the numerical solution to the solution of the continuous TDSE. Roughly speaking, stability has to be 
understood in the sense that the numerical solution will remain bounded. Or more precisely, in £^-norm and 
for all n, the stability condition is 

AxJ2m'=:\\rr;?\\r'-T- (201) 

j 

The finite difference scheme is also appropriate to solve TDSE in the velocity gauge. Indeed the transport 
operator can be discretized easily in a stable way as described below. First, it is recalled that the TDSE is 
given by 

dt^ = z9,,0 - tVc{x)i; + iA{t)d^i^. (202) 
The transport operator (hyperbolic) necessitates to upwind the discrete the operator 



A{tn)d,ijUdx,n ^ A{tn) ' /''\ ifA(t„)>0, 

ax 

Aitr,)d,i,{jdx,n ~ Aitn)-^^^-—^, ifA(t„)<0 

ax 



(203) 



The upwinding ensures the stability of the numerical scheme. The approximation of A(t)dx should then 
be done accordingly to the sign of A(i), or equivalently the approximation of the derivative in jdx is done 
accordingly to where the information comes from: from the left {dxtp ~ (V'J' ~ 4'^~i)/dx) if yl(t") > 0, 
and from the right [dxi^ ~ (V'j+i ~ i^j')/dx) if A{t'^) < 0. If this rule is not satisfied, the scheme becomes 
unstable, the solution blows up and as a consequence does not converge. 
In the length gauge, an order 2 (in space and time) scheme, on a point grid writes: 



(204) 

An I An+1 

The scheme can then be rewritten in the form 

^n+l^n+l ^ dt{B"ijj'' + F" + F"+i) , (205) 

where ^4"+^ and are sparse N x N matrices. Thus, the equation has become a linear system of equation. 
The latter can be stored on a computer by compressed sparse row storage [102) in order to avoid or limit 
storage issues. For sparse symmetric matrices A, the linear system can be solved by the conjugate gradient 
method, which is among the most efficient solvers. GMRES and Bi-conjugate gradient techniques are most 
efficient in non-symmetric cases [102) , which occur for instance when non- uniform spatial discretizations are 
used. The multi-dimensional TDSE can be solved similarly on a iV'^ point grids. As a consequence, for 
N large, this necessitates the storage of huge matrices and the solving of sparse linear systems. In order 
to limit the storage and computational time issues. Alternate Direction Implicit (ADI) methods are often 
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used, which mainly consists of sphtting the equation in each spatial direction and necessitate the solving 
of one-dimensional TDSE's. The consequence of this splitting is that the time step has to be reduced in 
comparison to direct methods to maintain a good accuracy. 



Finite volume method (FVM) for TDSE in the velocity gauge. We roughly describe how to derive 
a finite volume scheme for the TDSE. The interest of such a method is multiple. First it allows to consider 

any geometrical domains, with non-uniform cells (or volumes) as for the finite element method (the mesh 
can be designed to have finer cells in regions where the solution have strong gradients or singularities) . Then, 
this approach, based on a weak formulation of the equation, allows a very simple upwinding of transport 
operators (more generally hyperbolic operators) ensuring (binder condition on the time step) the stability 
(no blow-up of the numerical solution) and then, the convergence of the numerical scheme. 
For the sake of notation simplicity, we suppose that the physical domain il is a polygon (in 2-d) or polyhe- 
dron (in 3-d) . Domain fl is decomposed in cells or volumes (typically triangles in 2-d and tetrahedra in 3-d) 
denoted by Ki, then O = Uj-fCj. The starting point is the TDSE given by 

i^{r,t) = -^Ai,{r,t) + Vc{r)^{v,t) + iA{t) ■V^{r,t). (206) 
Then, integrating over the volume Ki, we get 

[ ilj(r,t)dT = — —[ Ailj{r,t)dr+ [ VJr)ip{r,t)dr + i [ div(A(t)Mr,t))dr{207) 

We denote i/'j^. the average of in Ki at time i", that is 

By the divergence theorem, where n^j is an outward normal vector to Ki and (Tj a surface on KiS boundary, 
we have 

iAt" At" 

= ^-^^2^^;^)^-^^^^^^'^^ ,,,,, 

-/^.F,(r)V(r,t)rfr, 



YO\{Ki 

where the time derivative was approximated by a finite difference. Because we are considering polyhedron 
volumes, the integration on the surface dKi can be rewritten as 

iAt" . , , 



At" , 



(210) 



where e^^^ denotes the faces (or edges in 2-d) of volume Ki. Now, V^" is an approximation of ijj^, and V^" 
an approximation of J^a-, A{t"')tjj{Y,t)dr/vol{Ki). Usual first order, explicit, finite volume schemes write 

iAf" 

= 1/." + -—Lf + At" J2 ) - i^t^VciVl^, (211) 

where Vc^i = Jj^, Vc{r)dr/vo\{Ki), and is an approximation to Atlj{r,t")dr/vol{Ki). Usually, L" is 
evaluated by reconstructing the Laplacian, using Ki^s neighboring cell values. The presence of Lf imposes a 
restrictive stability conditions on At" (typically in Atn ^ C( min^ vol(_ft')) ). The most important point to 
consider, and this is why FVS are well adapted to transport problem is the simple and stable approximation 
of the solution at the faces (or edges) of Ki. It is based on an upwinding process, which again ensures the 
stability of the scheme. More precisely, at a face e^'^ if we denote by Kj the cell sharing the edge e^*^ with 
Ki, then a stable approximation of , is given by 

_ r F," if A".n,„ >0, 
Kii) - j V^n jf A" • n,(.) < ^^^^^ 
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where n^(i) is the outward normal vector to the edge e^'^ of Ki. For instance, a naive approximation such 
as F") = {Vj"' + Vp)/2 would lead to numerical instability (for the same reasons as (12031) ). then to non- 
convergence. This is a natural extension of the approximation (|203p for FDM. To the best of our knowledge, 
FVM is not commonly used for approximating TDSE, although this approach has very nice computational 
and mathematical properties. In the length gauge, the use of FVM is less natural, but still applicable, due 
to the absence of transport term. Theory about FVM is very well developed, in particular for application 
in fluid dynamics and to a certain extent to electromagnetism. In the framework of quantum mechanics, 
however a lot has still to be done. 

6.2.4. Summary of numerical methods As seen in previous sections, the choice of numerical methods is 
closely related to the mathematical structure of the equation we plan to approximate. We have recalled 
above that transport operator approximation necessitates for stability reasons, an upwinding of the spatial 
derivatives. As a consequence finite difference or finite volume methods are perfectly adapted to this operator 
(Galerkin methods are by default centred techniques). However, the kinetic operator which is transformed 
into a symmetric bilinear form, by variational computations and is then well adapted to Galerkin approaches 
(finite element, spectral methods). To summarize the numerical solvers should be chosen accordingly to the 
fact that: 

• Kinetic operators: Galerkin methods (finite element, spectral methods). 

• Transport operator: upwind finite difference or volume methods. Galerkin's methods can also be used 
to approximate transport problem (using the Petrov- Galerkin Streamline Upwind, SUPG method). It 
consists of adding artificial numerical viscosity to stabilize the scheme. 

Note that in Coulomb gauge V • A = 0, it is possible to use a Lagrangian approach because in that case 

A • = div(A7/>) = 0. (213) 

We denote by {vi)i the set of grid points. The principle of Lagrangian methods consists of considering grid 
points as particles of a fluid propagating at velocity A. The kinetic equation 

idti}^-Mj + Vc{Y)^, (214) 

is solved by discretizing on the moving grid. That is at time we search for V'(r", ^n+i), that is the solution 
at time in+i defined on the grid points located in r" at time i„ and moving at velocity A". This technique 
is in particular appropriate to the acceleration gauge. 

All the presented methods can usually be coupled with mesh adaption. Mesh adaptation which can be 
based on wavelet decomposition [103] or local error estimators [104| , AMR |105) , is a very useful tool from a 
practical point of view. These are technical methods that dynamically adapt the spatial mesh, in function 
of the solution regularity: where the solution have large gradients or singularities, degrees of freedom are 
dynamically added; where the solution has slow spatial variation, degree of freedom are removed. In fine, 
mesh adaptation allows for very accurate numerical solutions with acceptable computational times and data 
storage. 

6.3. Boundary Conditions 

The boundary condition question is essential in the discretization of TDSE's. Due to the intensity of the 
laser pulse, the wave function is extensively delocalized necessitating a large computational domain. In 
multidimension and for multi-electron systems it is important to choose appropriate boundary conditions 
to avoid spurious reflection of the wave function at the computational domain boundaries. From a gauge 
invariance point of view, it is also crucial, as non-exact boundary conditions will necessarily lead to a 
discrepancy of the gauge invariance. As far as we know, no study exists on the effect of (the choice of) the 
boundary conditions on the gauge invariance. Taking Dirichlet or Neumann boundary conditions leads to 
important numerical oscillations and reflections at the boundary of the domain, interacting with "physical" 
waves inside the domain. Even if this kind of methods allows effectively to reduce spurious reflections, there 
are often empirical (see for instance [106] in this framework), as some "parameters" have to be adapted for 
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each numerical situation. Outside the bounded domain, the Coulomb potential is assumed to be negligible 
and the laser-molecule TDSE can then be solved "exactly" using for instance the Volkov state propagator 
(see [106j ). Ideally we would like to impose boundary conditions such that the solution in the whole space 
restricted to a bounded domain is equal to the solution in fl (that is without spurious reflections). Then 
outside this fictitious domain the wave function is accurately approximated and can be updated using for 
instance Volkov state propagator, or by another TDSE solver associated to other nuclei. As an illustration 
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Figure 3. Electronic wave function delocalization 



of this issue let us consider the following case: Let us consider the following simplified model without electric 
field: 

1 

2' 



idtu{y, t) + ^d^yu{y, t) - K(y) • u(y, t) = 0, 



u(y,0) 'uo(y) 

One considers the domain 51 x [0,T] and denotes by F the boundary of Vt. One then looks for v solution of 

idtv{y,t) + ]flyv{y,t) - V^iy) ■ v{y,t) = 0, y G f2, 

B{y,dy,dt)v{y,t)=0, y e F, 

I v{y,0) = UQ{y), y er 
such that 



(215) 

The main problem consists then of finding an adequate (pseudo-)differential boundary operator S on F such 
that pi5|) occurs, see Fig. ID As is well known these conditions, called Neumann-Dirichlet (or Dirichlet- 
Neumann), see Appendix A, are non-local in time (and in space in higher dimension). Denoting by n the 
outward normal of F and dn is the trace operator on F we obtain: 



idMy.t) + -dlv{y,t) - VM ■ v{y,t) = 0, y G 11, 



v{y,t) = -eW4V2/,* ^""/.''"\dr, y e F. 



- r) 
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U| = V 



Figure 4. Transparent boundary conditions 



This approach has been very well described in particular in |107) . and some results can be found in [108) . 
1 109) or [llOj . We also refer to for the first presented discretization of non-local transparent boundary 
conditions for TDSE's. As unfortunately these conditions are non-local in time, many attempts have been 
devoted to find efficient numerical approximations of these conditions. 

To illustrate this technique, we propose again a simple benchmark. We suppose that the Coulomb 
potential is equal to zero. 

r 1 

idtv{y,t) + -dyyv{y,t)~yE{t)v{y,t)=Q, ye [-10,10], i ^ 0, 

^ v{y,Q) = vo{y) = e^'ye-y\ 

The benchmark we propose is as follows. The fictitious domain is f2=[— 5,5]. We impose the Neumann- 
Dirichlet boundary conditions coupled with the laser as described above, at x-r — —5 and a;r — 5. We 
compare our numerical solution with the solution obtained using Dirichlet boundary conditions at a;_r = ^5 
and xr = 5 and with a reference solution obtained numerically on a large domain. 

The "Neumann-Dirichlet numerical (fig. [5]) solution" is then far less reflected (even if a little reflection exists) 
than the "Dirichlet solution" . Here, note that the grid is coarse and small, so that the influence of spurious 
reflections can be obviously diminished using a larger grid and smaller space steps. We also represent the 
norm error between the reference solution and the Dirichlet and Neumann-Dirichlet solutions. 



7. Gauge transformation and the Dirac equation 

In the last few decades, laser intensities achieved in experiments have increased due to new technical advances. 
It is now possible to consider laser fields with / ^ lO^'^W/cm^, and higher |112| . In this new regime, the 
electron starts to move at relativistic velocities. For instance it is demonstrated in '113| , by using an argument 
based on classical relativistic mechanics, that when 

^ > 1 (216) 

where Eq is the electric field and a; is the laser frequency, an electron is accelerated to a regime where 
relativistic effects start to be important. The mathematical description of an electron subjected to such 
intense electromagnetic fields necessitates a relativistic treatment |114| 1113) and thus, theoretical efforts 
should be based on the Dirac equation instead of the non-relativistic Schrodinger equation. 

Most of the discussion concerning gauge invariance can be applied to the Dirac equation, which gives a 
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- Neumann-Dirichlet 




Figure 5. Comparison between the reference solution and the numerical solutions obtained with Dirichlet 
and Neumann-Dirichlct boundary conditions 



quantum relativistic description of the electron. The latter is given by jll5) 

idtip{t,x.) = |a • [p- A(t,x)] + (3mc^ + el4J7(i, x)|v'(t, x), (217) 

where ipitjX.) G L^(M^) (8) is the time and coordinate dependent four-spinor, I„ is the n by n unit matrix 
and ai,/3 are the Dirac matrices. This equation describes physically the relativistic dynamics of a single 
electron (spin-1/2) subjected to an external electromagnetic field. As in the Schrodinger equation, the latter 
is introduced by using the minimal coupling prescription which allows to preserve the gauge invariance of 
the equation: this will be shown below. In this work, the Dirac representation is used where 

"0 

Cti '■ = 







I2 
-I2 



The cr,: are the usual 2x2 Pauli matrices defined as 

and cr. 



(218) 








1' 




-i 




1 





, <7y := 


i 



-1 



(219) 



There exist many other representation of Dirac matrices as they are defined abstractly by their anti- 
commutation relations. A list of representation can be found in |116| I115j . They are related to each 
other by unitary transformations. 

To show the gauge invariance of the Dirac equation, it is possible use the Lagrangian formulation and 
then, to demonstrate that the Dirac Lagrangian obeys the symmetry condition under a gauge transforma- 
tion. Here, we will use the equation of motion. We have two arbitrary gauges where the wave function and 
gauge field are related by 

V-(i)=e'V\ (220) 
A(i) ^ A(2) + VF, (221) 
[/(I) = _ dtF. (222) 

The wave equation in gauge 1 obey the following Dirac equation: 

(223) 



(224) 



idtij^^\t,x) = ^a- ^p- A(i)(t,x)J + (3mc^ + el4U^^\t,yi)ji;^^\t,x.) 
Using the gauge transformation to gauge 2, we get 

idti^'-^'^it,x) = \ a- [p- A(2)(t,x)j + /3mc^ + eIiU^^\t,y:)\ij^^\t,x) 
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and thus, the Dirac equation is invariant under gauge transformations. This is the same resuh as for the 
Schrodinger equation considered previously. Therefore, many of the results presented earlier can also be 
applied to the Dirac equation. 
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Appendix A. Exact solutions 

In many situations (more or less physical) the laser-molecule TDSE can be solved explicitly or approximately. 
We here shortly recall some of these particular situations. 

• Volkov: Exact solution to potential-free (vacuum) TDSE or TDDE can be obtained by solving 
analytically. 

idti^^^{p-A)^i: (A.l) 
Zm 

The principle is as follows. We apply a Fourier transform in space, which leads to a differential equation 
which is can be easily solved. Taking the inverse Fourier transform leads to the Volkov wavefunction: 

■i/'(r,t) = ^^exp ir^ E(s)ds - {J 'E{T)dTfds 

fn(fnE(r)dT)ds- r')% , , , , 

V'o(r')exp(^ ^"^-^^ ' ^jd'r', (A.2) 

where E denotes the electric field. From the numerical point of view the main issue comes from the 
non-locality of this solution. Naturally this method is presented in any quantum physics book and can 
be extended to relativistic situations (Dirac equation). 

• Dirichlet-Neumann. Laser-free TDSE in vacuum can be solved analytically in length gauge. We 
shortly recall the principle here. From the laser-free Schrodinger in vacuum (in 1-d) 

Tpt- idxx'il' ^ {dV^ - e'fda:){dl''^ + €"^3^)11} ^Q, V^a^, i) ^I^Koo (A.3) 

involving the Dirichlet-to-Neumann pseudo-differential operator, we can solve explicitely this equation 
(under the Somerfeld radiation condition) 



yjTT{t - t) 

This is naturally a fundamental result to derive analytical solution to TDSEs. It is also a very used tool 
to derive absorbing boundary conditions [1171 11181 1119) . 

• Lewenstein Model/SFA. This approach (which was already discussed above), although not exact, 
allows to find solutions to TDSE for intense laser pulses, including the Hamiltonian continuum |12) . 
That is we search for wavefunctions of the form: 

i/'(r, t) — exp {iIpt)(a{t)^o + j d^pa(p,t)i^p) (A.4) 

which corresponds to a decomposition of the wavefunction on the ground ■00 E^nd continuous ^-p states. 
Ip denotes the ionization potential, the density h is the amplitude of the continuum state ip-p. This 
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model has been extensively developed and validated as it gives a very good description of multiphoton 
ionization and high order harmonics generation. The dipole can be calculated as follows 

d{t) = J^,\i;{r,t)\^rd^r 

(A.5) 

= i£J dt'd^pE cos(t')rfx(p - Mt')d* (p - A{t)) exp ( - ^^(p, t, t')) + ex. 

where 

S{p,t,t') = j'dt"(^^^^—^^^ + Ij,). (A.6) 
is the action of the free electron in the vector potential A{t) 
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